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ABSTRACT 

Computer simulation techniques were used to determine 
equilibrium positions and binding energies of inert gas atoms 
implanted in a tungsten crystal and to investigate the po- 
tential wells around these equilibrium positions in both per- 
fect lattices and relaxed lattices. Stable positions were 
found for inert gas interstitials near lattice atoms in the 
third and fourth layers of the crystal. Interstitial posi- 
tions near atoms in the first and second layers of the crys- 
tal appeared to be unstable if they exist at all. As a 
result of potential well studies, it was concluded that the 
mechanism associated with equilibrium position formation was 
a combination of local liquefaction of the lattice structure 
and interaction of the interstitial with lattice atoms. 
Equilibrium positions were found to be ill-defined regions 
in the general (110) direction. The binding energy deter- 
mined for an interstitial site near a lattice atom in the 
third layer of the crystal was in excellent agreement with 
experimental results. 



2 



TABLE OF CONTENTS 



I . INTRODUCTION 7 

II. THE NATURE OF THE PROBLEM 9 

A. THE LATTICE DYNAMICS PROBLEM IN THE COMPUTER 9 

B. HISTORICAL DEVELOPMENT OF COMPUTER SIMULATION 

OF LATTICE DYNAMICS 11 

III. THE SIMULATION PROCEDURE 16 

A. ' THE MODEL 16 

1. The Crystal 16 

2. The Potential Function 19 

3. The Timestep 21 

a. General Discussion of the Timestep 21 

b. The Average Force Method and Timestep 

Determination 23 

c. Procedures Used to Determine DTI 24 

B. THE PROGRAMS 26 

1. The Static Program 26 

2. The Dynamic Program 28 

3. Equilibrium Sites Viewed As Potential 

Wells 29 

4. Determination of Possible Interstitial 

Sites 30 

a. Implantation Procedure 30 



3 



IV. PRESENTATION OF DATA 31 

A. STATIC SIMULATION 31 

1. Initial Simulations 31 

2. Site 39 31 

3. Site 14C 34 

4. Force Calculations 35 

5. Investigation of Potential Wells 36 

a. Potential Well Studies in a Perfect 

Lattice 36 

(1) Interpretation of Results 39 

b. Potential Well Studies in a Relaxed 

Crystal 40 

(1) Interpretation of Results 41 

B. DYNAMIC SIMULATIONS 43 

V. CONCLUSIONS AND RECOMMENDATIONS 45 

APPENDIX A: SUBROUTINE PLUCK 47 

APPENDIX B : COMPUTER PROGRAM GLOSSARY 50 

APPENDIX C: ERROR IN THE LITERATURE 59 

COMPUTER PROGRAMS 65 

BIBLIOGRAPHY 97 

INITIAL DISTRIBUTION LIST 99 

FORM DD 147 3 100 



4 



LIST OF TABLES 



I. Interstitial and Lattice Atom Displacements 

From Reference Site 89 32 

II. Results of Replacement Runs in Site 89 38 

III. Interstitial Potential Energies for Site 89C 38 



5 



LIST OF DRAWINGS 



1. Atom Numbering Procedure 61 

2. The Function of Subroutine Pluck 62 

3. Xenon Behavior in a Perfect Lattice 63 

4. Argon Behavior in a Relaxed Lattice 64 



6 



I. 



INTRODUCTION 



In 1968 Kornelsen and Sinha [1] of the National Research 
Council of Canada published results of radiation-damage ex- 
periments performed on tungsten. In these experiments a 
tungsten surface was bombarded with ions of neon, argon, 
krypton, and xenon respectively. Then, while the tungsten 
was heated, gas desorption rates were measured as the gas 
evolved from the crystal. The resultant desorption spectrum 
was interpreted to yield a binding energy spectrum for the 
trapped particles. In 1970 Professor Don E. Harrison, Jr., 
of the Naval Postgraduate School undertook the modeling of 
these experiments utilizing computer simulation techniques 
in order to provide a means for interpretation of their re- 
sults. Two successive thesis research efforts [2,3] have 
been specifically directed toward the investigation of inert 
gas implantation in a tungsten crystal. It was anticipated 
that a corollary to the successful computer simulation of 
this problem might possibly be an improved understanding of 
the interstitial atom stabilization mechanism in tungsten, 
with more general application to other materials. 

The investigations reported in this paper are a con- 
tinuation of work begun by Vine [21 and Tankovich [3] under 
Harrison's supervision. The simulation procedures followed 
were a combination of static and dynamic approaches to the 
problem. The static portion of the problem entailed 
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interstitial implantation of an inert gas atom in a tung- 
sten crystal and the subsequent relaxation of the crystal 
until the equilibrium position of the interstitial could be 
ascertained. In the dynamic portion of the problem de- 
creasing amounts of energy were imparted to the interstitial 
in its equilibrium . position until the minimum energy, and 
direction, which still allowed the interstitial to escape 
from the crystal was determined (i.e., the binding energy 
of the particle) . The binding energies thus determined 
provided a basis for comparison with the results of Kornelsen 
and Sinha. 
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II. THE NATURE OF THE PROBLEM 



A. THE LATTICE DYNAMICS PROBLEM IN THE COMPUTER 

For a little over a decade high speed digital computers 
have exhibited their usefulness as tools for investigation 
of physical systems. Specifically, the inherent periodicity 
and order of crystal lattices have made the study of lattice 
dynamics particularly adaptable to investigation by computer 
simulation. It is a relatively simple matter to "construct" 
a crystal of the desired body-centered cubic or face-centered 
cubic structure in the computer. Various types of point de- 
fects can also be "created" in the crystal with relative 
ease. A vacancy is obtained by simply removing an atom from 
the crystal, while interstitials can be created by implanting 
an additional atom in the crystal. Two types of interstitials 
have been used in investigations of lattice dynamics, atoms 
identical to the crystal lattice atoms ("self-interstitials") 
and atoms different from the crystal lattice atoms 
(interstitials) . 

Although a crystal lattice containing a point defect can 
be easily represented in a computer, modeling of the dynamics, 
which allows alterations _of the crystal structure resulting 
from the presence of the point defect, is more involved. 

The dynamic portion of the problem of computer modeling of 
lattices can be characterized by four hey decisions which 
must be made. 
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•First, it is necessary to find some mathematical re- 
lationship to govern the interaction among the atoms of the 
crystal. This interaction is a complex many body problem 
which is nearly always approximated in computer simulations 
as a sum of appropriate two body interactions. With this 
approximation it is then possible to represent these two 
body interactions by some type of potential function, which, 
in turn, can be used to determine forces on individual atoms. 

Secondly, since the number of calculations required is 
directly related to the number of atoms in the crystal, a 
judicious choice of crystal size must be made. Large crystals 
would give more accurate results, but would also require more 
computational time. 

A third key problem which must be solved results from 
the inability of a digital computer to perform a direct in- 
tegration. Since the lattice dynamics problem is most often 
directed at a determination of lattice atom positions after 
some type of interaction, an integration of the equations of 
motion is required. A choice of the numerical integration 
technique to be used must therefore be made. 

Last, and probably most important, since an iterative 
process is used to integrate the equations of motion, the 
length of the time interval of the iteration must be chosen 
so that the force variations within this time interval are 
small and, consequently, stability of the system is insured. 

At first glance, this suggests that only very small timesteps 
should be taken, but this would require excessive computer 
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time. Consequently, some variational method of timestep 
determination would be optimum. 

B. HISTORICAL DEVELOPMENT OF COMPUTER SIMULATION OF 

LATTICE DYNAMICS 

The four basic problems of computer simulation of lattice 
dynamics have been solved in various ways in previous 
investigations . 

In 1960 Gibson, Goland, Milgram, and Vineyard [4] 

(referred to hereafter as Gibson, et.al) of the Brookhaven 
National Laboratory published the first complete statement 
of computer simulation procedure as a tool for investigation 
of radiation damage of crystal materials. Due to the com- 
plexity of the radiation damage problem and the inadequacy 
of analytical methods as a means of analysis of the damage 
processes, Gibson, et.al. turned to a numerical integration 
of the problem utilizing high speed computers. Their pio- 
neering work gives insight into such pertinent subjects as 
crystal size, choice of potential functions, computational 
methods for solving the equations of motion, and timestep 
determination . 

Specifically, their research utilized copper as the 
crystal material since its relatively simple face-centered 
cubic structure and its widespread use in actual experimen- 
tation made it particularly appropriate for initial investi- 
gation. A Born-Mayer repulsive potential was used to describe 
atom-atom interaction, and a constant cohesive force was ap- 
plied to each atom on the crystal surface to balance the 
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Born-Mayer repulsion. A central difference iterative 
procedure was used to integrate the equations of motion. 

The procedure used for timestep determination is of parti- 
cular significance since the principles involved are still 
the governing criteria for choice of timestep duration. The 
fundamental observation was made that the greatest stress 
on the crystal system was a result of the strongest atom- 
atom interaction in the system. This interaction was there- 
fore chosen as the basis for timestep determination. Simply, 
the timestep duration was chosen to be inversly proportional 
to the velocity created by the strongest atom-atom 
interaction. 

In addition to verifying the applicability of computer 
simulation techniques to radiation damage studies and pro- 
viding specific information on collision chains and focusing 
phenomena in crystallites, the work of Gibson, et.al.[4] 
determined that the only stable configuration for self- 
interstitials in a face-centered cubic crystal was the <100 ) 
split configuration. The split configuration implies that 
the interstitial causes a lattice atom to split away from 
its normal lattice position, and then both atoms, the inter- 
stitial and the lattice atom, share (or split) the normal 
lattice site when equilibrium is reached. 

Johnson and Brown [5] confirmed the split configuration 
as the only stable interstitial position in their studies 
of copper utilizing a Born-Mayer repulsive force between 
nearest neighbors and an elastic continuum containing the 
remainder of the atoms of the crystal. Erginsoy, Vineyard, 
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and Englert [6] (referred to hereafter as Erginsoy, et.al.) 
extended these calculations to the body-centered cubic case 
by performing investigations using « iron. Computational 
techniques paralleled those of the Brookhaven group except 
for the choice of a potential function. After experimen- 
tation with a Born -Mayer potential and a Morse potential 
with parameters derived by Girifalco and Weizer [7] , 

Erginsoy, et.al. [6] settled on a combination potential which 
utilized an exponentially screened Coulomb potential for 
close approaches, a Born-Mayer potential for first nearest 
neighbor, and a Morse or modified Morse potential for 
higher order neighbors. This research verified the split 
configuration as the only stable configuration for body- 
centered cubic structures, but the orientation of the split 
was found to be in a (110) plane, (vice a (100) plane as in 
the face-centered cubic case ) . R.A. Johnson [8] confirmed 
the split interstitial orientation in similar work with cc 
iron, vanadium and tungsten. 

In their work on collisions between a copper atom and 
a copper lattice. Gay and Harrison [9] introduced the average 
force procedure for integrating the equations of motion. A 
complete description of this procedure has been reported by 
Gay, Effron, and Harrison [10]. In two more recent efforts 
Johnson and Wilson [11,12] determined new potential functions 
for use in face-centered cubic and body-centered cubic defect 
calculations and published results of defect calculations of 
helium in various metals. 
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Parallel to the computer simulation efforts described 
above, Kornelsen [13,14] and Kornelsen and Sinha [1,15] 
have published the results of extensive experimentation 
involving the interaction of inert gas ions with tungsten. 
Under Harrison's supervision. Vine [2] initiated attempts 
to devise a computer model of Kornelsen 's experiments on 
neon implantation in tungsten. He attempted to develop a 
relationship between equilibrium positions of tungsten 
interstitials in a tungsten crystal determined, first, by 
using a Born-Mayer repulsive potential to describe tungsten 
interstitial interaction with a tungsten lattice then, by 
using the same tungsten-tungsten composite potential which 
was assumed between atomsof the tungsten crystal. This re- 
lationship would then have been applied to neon equilibrium 
positions derived from a repulsive potential to provide a 
comparison with Kornelsen' s data. These efforts met with 
little success. 

Follow-on work by Tankovich [3] utilized the static/ 
dynamic approach described in more detail in Section III-B. 
The potential function used in these investigations was a 
composite Born-Mayer and Morse potential joined by the best 
cubic fit in the area of intersection which had previously 
been developed by Harrison and Moore [16] . Static program 
runs confirmed the (110) split interstitial for helium, 
neon, argon, krypton, and xenon point defects in a tungsten 
crystal for possible interstitial sites in planes three 
through six of the ten plane crystal used. Preliminary 
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dynamic testing of an argon equilibrium position in plane 
four of the crystal was begun with limited success. 

Of particular significance in Tankovich's work [3] was 
the introduction of a timestep decrementing process into the 
static program. In previous investigations the timestep 
duration was chosen at the outset of the problem and remained 
constant through all computations. At this author's sug- 
gestion, a timestep decrementing process was devised which 
allowed a more rapid approach to the equilibrium positions 
with a concomitant saving of computer time required for 
computation. (See Section III-A-3 for a discussion of the 
timestep.) 
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III. THE SIMULATION PROCEDURE 



A. THE MODEL 

1 . The Crystal 

Two factors are of primary importance in determining 
the size of the crystal to be used in computer simulation 
investigations. First, the crystal must be large enough to 
provide realistic results, and, second, the crystal must be 
as small as possible in order to minimize computer time re- 
quired for calculation of atom- atom interactions. 

After experimenting with crystals of various sizes, 
Tankovich [3] determined that a crystal with ten planes of 
atoms in each coordinate direction, which is equivalent to 
five unit cells, (10 x 10 x 10) was suitable for static in- 
vestigations with tungsten. The same crystal was consequently 
used in the investigations reported in this paper. The lat- 
tice unit, or distance between adjacent (100) planes of 
atoms, for a body-centered cubic tungsten crystal is 1.58$L 
All distances in the computer simulation were measured in 
lattice units. The lattice constant for the tungsten crystal 
is 3. 16& (or two lattice units) , and the nearest neighbor 
distance is \| 3 lattice units. 

The same numbering sequence for atoms of the crystal 
that was employed by Tankovich [3] was used in these investi- 
gation. (See Figure 1.) Atom number one was always assigned 
to the interstitial atom. A rectangular coordinate system 
was placed on the upper, left hand, front face of the crystal. 
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The positive "x" direction was chosen to the right from the 
origin, the positive "y" direction was down, and the positive 
"z" direction was to the rear. Atoms were numbered consecu- 
tively, beginning with atom number 2 at the origin and con- 
tinuing until all atoms of the Y = 0 plane had been numbered. 
This same numbering sequence was then followed for each Y 
plane of the crystal until all 250 atoms of the crystal had 
been numbered. Figure 1 shows the numbering of the atoms 
in the Y = 0 and Y = 1 planes. The numbering of atoms in the 
remainder of the planes follows the same pattern. 

For the dynamic program it was felt that a smaller 
crystal could be used and still yield meaningful results. 

The rationale for this determination was based on the go-no 
go (i.e., escape or no escape) character of the dynamic pro- 
gram. The dynamic program essentially provides an answer to 
this question - Will an interstitial atom escape from the 
crystal if it is given a specific energy and directed in a 
specific direction? If the energy dissipation mechanism 
provided by collisions with, and close approaches to, the 
lattice atoms along the path traveled is great enough, the 
interstitial will remain in the crystal. If these energy 
dissipation mechanisms are not large enough to overcome the 
kinetic energy of the interstitial, the interstitial will 
escape. Since the atoms of the crystal along this line of 
motion must provide the mechanism to dissipate the inter- 
stitial ' s kinetic energy, for escape to be prevented, only 
atoms in the horizontal planes above the interstitial and 
atoms in vertical planes within a few lattice units of the 
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interstitial will affect the possible escape. The remainder 
of the atoms of the crystal would not have time to react 
with the interstitial or affect its movement. As a minimum, 
the smaller crystal could be used to eliminate excessively 
high or low energies from consideration at a considerable 
savings of computer time and, instead, give a limited range 
of energies to be checked by dynamic runs using the entire 
crystal . 

To implement this procedure SUBROUTINE PLUCK was 
developed. (See Appendix A for a complete discussion of 
SUBROUTINE PLUCK.) Basically, SUBROUTINE PLUCK uses the re- 
sults of the static program, but causes a crystal to be 
printed out that only contains the interstitial and all lat- 
tice atoms from two planes below the interstitial to the 
surface plane of the crystal and all atoms in vertical planes 
which are within two (or three) lattice units of the vertical 
plane containing the lattice site shared by the split inter- 
stitial. (See Figure 2.) The savings in computer time re- 
sultingfrom the use of this smaller crystal is demonstrated 
by the following comparison. A thirty timestep dynamic run 
with the entire 10 x 10 x 10 crystal (250 atoms) used approxi 
mately three and one half minutes of computer time. A thirty 
timestep dynamic run using a 7 x 5 x 7 crystal (60 atoms) 
used slightly less than one minute of computer time - a 71% 
savings in computer time I 

Most of the dynamic runs made' during the course of 
these investigations utilized the 7x5x7 "PLUCK" crystal. 

Final confirmation of minimum energies was determined using 
the entire crystal. 
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2. The Potential Function 



As mentioned previously, the many-body interaction 
which characterizes actual lattice dynamics is approximated 
in computers by many two-body interactions. These two-body 
interactions are represented in the computer by some type of 
central, pairwise potential function. Various types and 
combinations of potential functions have been considered for 
use in computer simulation investigations of lattice dynamics. 

The choice of the potential function must be made 
with due consideration to the range of applicability of the 
potential function, the correlation of potential function 
parameters with observable properties of the material being 
investigated, and the amount of computer time required for 
calculations using that potential function. The potential 
function used in these investigations was the composite 
Born-Mayer potential and Girifalco and Weizer Morse potential 
used by Tankovich [3,16] . This composite potential is con- 
structed as follows: 

a. Region 1- (r < 1 . 5&) 

The atom-atom interaction at close approach is 
represented by a Born Mayer repulsive potential of the form, 

= exp (A+B r i j) (1) 

where is the interaction energy between particles i and 

j and r^j is the distance between particles i and j . 

b. Region 2-(1.5& < r < 2.0&) 

This portion of the potential function is ob- 
tained by computing the best cubic equation between the value 
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of the Born Mayer potential at 1.5& and the value of the 
Morse Potential at 2.0&. 

c. Region 3-(2.oX < r < 5.38&) 

For equilibrium and greater separations, a 

Morse potential of the form, 

$ ij = D [ ex P£ _2 a (r ij -r 0 ) }“ 2 exp{- a (r i j-r Q ) }] (2) 

where is the interaction energy between particles i and 

j, D is the dissociation energy of the particles i and j, r^j 
is the distance between particles i and j, and r Q is the 
equilibrium separation, is used. The Morse potential was 
computed so that the tail of the function was truncated to 
zero at 5.38&. This effectively meant that atoms out to the 
fourth nearest neighbor were included in calculations of 
interaction energies. Girifalco and Weizer [7] had pre- 
viously determined Morse potential parameters for tungsten 
and other elements which included interactions out to the 
150th nearest neighbor. Use of the complete function, 
however, would have required an excessive amount of computer 
time for calculation. Additionally, contributions to the 
interaction energy of all atoms beyond the fourth nearest 
neighbor is essentially insignificant for our calculations. 
Many computer simulations of lattice dynamics of body- 
centered cubic materials have utilized potential functions 
which only included first and second nearest neighbors with 
satisfactory results [6,11]. To check the adequacy of this 
potential function in describing this crystal system, the 
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largest binding energy observed in the crystal model 
(-8.283 eV) may be compared with the experimentally de- 
termined heat of sublimation for tungsten (-8.8 eV) reported 
by Harrison and Magnuson [17]. This agreement within 5.9 
percent was considered satisfactory for these simulations. 

3 . The Timestep 

a. General Discussion of the Timestep 

Most of the early simulations of lattice dynamics 
used a central difference method as the numerical procedure 
of integrating the equations of motion of the atoms in the 
crystal. See Gibson, et.al [4] or Gay, Effron and Harrison [10] 
for an explanation of the central difference method of numer- 
ical integration. The investigations reported in this paper, 
however, used the average force method which is completely 
described by Gay, Effron, and Harrison [10] . Inherent to 
both numerical methods of integration is the replacement of 
time derivatives in the equations of motion with a finite 
time difference, i.e., the timestep interval. As mentioned 
previously, the strongest interatomic interaction places the 
greatest demand on the system; consequently, the timestep 
duration is usually determined through consideration of the 
strongest interatomic interaction of the system. These in- 
vestigations followed the procedure of Gay, Effron, and 
Harrison [10] and determined the timestep duration by a con- 
sideration of the maximum displacement that the most ener- 
getic atom of the crystal was allowed to move. This parameter, 
referred to hereafter and in all computer programs as DTI, is 
determined as follows; 
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atom of the 



(1) The equation of motion of i^ 
crystal during timestep interval AT can he written in the 
form 

x i (t+AT)= x i (t) + [v i (t)+ <F i >AT/2m]AT (3) 

or in the equivalent form 

AXf = (v ± + <F i >AT/2m) AT. (4) 

(2) Rearranging terms of equation (4) yields 

AT = Ax ± /(v i + <F i >AT/2m) (5) 

where AT is the timestep interval, Ax^ is the displacement 
tTi 

of the i atom during the timestep, v^ is the velocity of 
the i^ 1 atom, and (F^) is the average force on the i^ 1 atom. 

From equation (5) it can be seen that the 
timestep interval is a function of both the kinetic energy 
and the force. Rather than solving this quadratic equation 
for the timestep interval, the average force method considers 
separately the cases when energies dominate forces 
(v^»(F^ ) AT/2m) and when forces dominate energies 
(<F^)AT/2m » v i ) . Solutions for each of these cases yields 
a different value for the (next) timestep interval. 

In energy dominant cases, 

AT = Ax^/v^ (6) 

or, expressed in the form of energies, 

AT = Ax i (m/2T m ) ^ (6a) 

and, in force dominant cases, 

AT = (2m Ax i /<F i >) J5 , (7) 



22 



where of equation (6A) represents the kinetic energy 
of the atom of the crystal with the greatest kinetic energy, 
and (F^) of equation (7) represents the average force on 
the atom with the maximum force, the Ax^ of each equation 
becomes DTI. 

In these investigations, DTI was set at the 
beginning of the program. Ideally, the proper choice of 
timestep duration for the first timestep and the proper 
choice of DTI would lead to a smooth movement of the inter- 
stitial to its equilibrium position in a manner similar to 
the movement of a critically damped oscillator. 

It was anticipated that energies would 
dominate forces in early timesteps which would lead to 
timestep determination from the energy equation, equation 
(6A) . At some point in the crystal relaxation procedure, 
energies should have been dissipated to the point where 
further timestep determination would become force dependent 
(equation (7) ) . 

b. The Average Force Method and Timestep Determination 
In the average force method of integration of the 
equations of motion, velocities (i.e., energies) of, and 
forces among, all atoms of the crystal are computed with the 
atoms in their initial positions. Based on these forces and 
energies and the initial timestep duration, new positions 
for all atoms of the crystal are computed. The forces at the 
new positions are then averaged with the forces at the ori- 
ginal positions to determine the average force, and hence the 
final positions of all atoms at the end of the first timestep. 
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In the meantime, the maximum force exerted on an 
atom of the crystal in its original position and the maximum 
force exerted on an atom of the crystal in its final (aver- 
aged) position are used with the DTI and equation (7) to 
calculate two possible alternatives for the next timestep 
interval. Likewise, the energies of the most energetic atoms 
in both original and final positions are used with DTI and 
equation (6A) to calculate two other possible alternatives 
for the next timestep interval. These four alternatives are 
then compared, and the smallest is chosen as the next time- 
step interval. 

c. Procedures Used to Determine DTI 

Vine[2] utilized a constant DTI in all of his 
computations. Since the movement of an interstitial to an 
equilibrium position cannot, in general, be characterized by 
a small range of energies and forces and since DTI should be 
closely correlated with energies and forces at least in ap- 
propriate regions, the use of a constant DTI in all calcu- 
lations made the initial choice of DTI extremely critical. 
Success could only be attained by resorting to small DTI 1 s 
and concomitant excessive program run times (s 100 times teps) . 

Tankovich [3] obtained more satisfactory results 
by successively decrementing DTI for each timestep. This 
procedure allowed long timesteps in early portions of a run 
when the interstitial was far from its equilibrium position. 

As equilibrium was approached, the decrementing process had 
progressed to the point such that significantly smaller and 
smaller timesteps were taken allowing a smooth arrival at the 
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equilibrium position. Additionally, this was accomplished 
at a considerable savings of computer time (~ 30 timesteps) . 

Although this decrementing process for DTI pro- 
vided considerable improvement, on occasion, the final few 
timesteps used such a small DTI that practically no atom 
movement was discernible. During these investigations this 
situation was alleviated by incorporating a minimum DTI into 
the decrementing process. This insured that atom movement 
was still discernible near equilibrium and assisted in 
guarding against possible false assumption of equilibrium 
because of the relatively small movement observed under the 
continuous decrementing process. 

In an attempt to more fully understand the 
mechanisms of the static solution, the computer program was 
adjusted so that the maximum force, the atom upon which this 
force was exerted, and the four "new" possible timesteps for 
each timestep interval were printed out after each run. It 
was observed that the timestep calculations based upon "new" 
and "old" energies were overwhelmingly the basis for timestep 
determination. To insure that the force dependence of the 
timestep was not being unduly disregarded, the program was 
adjusted so that DTI was determined solely as a function of' 
the minimum of the two forces. These "force calculations" 
of equilibrium positions agreed with "energy calculations" 
within 0.03 lattice unit. 
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B . THE PROGRAMS 



Although the basic computational procedures contained 
in the computer programs for the static and dynamic portions 
of the problem were essentially the same, an understanding 
of the subtle differences between the programs and an ap- 
preciation for several computational tools and procedures 
used in the programs have to be gained prior to further dis- 
cussion of the actual investigation and results. 

1 . The Static Program 

In the static program, a tungsten crystal of 
appropriate size was created in the computer. An inter- 
stitial atom was then implanted at a chosen site within 
the crystal. Potential energies and mutual forces of all 
atoms were computed. The crystal was then allowed to relax 
in appropriately chosen timesteps. At the end of each time- 
step an energy dissipation mechanism was introduced in the 
form of a predetermined damping factor which was used to 
decrease all velocity components of the atoms of the crystal. 
The next timestep interval was then computed, and the process 
was repeated until an indicated number of timesteps had been 
completed. If an equilibrium position had been reached, 
positions and energies of all atoms in the crystal, including 
the interstitial, and other pertinent data could be punched 
out on cards for later use in the dynamic program. 

a. Equilibrium Positions of Interstitials 

In interpreting the results of the static simu- 
lations it was necessary to arrive at some criterion to use 
as a determining factor for the interstitial's final 
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equilibrium position. It was expected that the crystal 
would relax around the split interstitial site forming a 
"pocket" within the crystal which would interrupt the 
periodicity of the lattice. Stability in this configuration 
was determined to have been reached when the atoms of the 
deformed crystal, including the interstitial, had been 
allowed to relax (i.e., adjust to the presence of the 
interstitial) to the point where their kinetic energies 
were all below thermal. (< 0.025 eV) 

' If more than one possible equilibrium position 
met this criterion, it was felt that different positions 
within the "potential well" of the equilibrium site were 
probably being observed. The position in which the inter- 
stitial atom had the smallest amount of potential energy was 
then chosen as the equilibrium position for the lattice site 
under investigation. 

b. Handling of Oscillations Near Equilibrium 

While performing the "force calculations", some 
runs began with the interstitial moving toward an apparent 
equilibrium until an oscillation, or "rattle", developed 
about the suspected equilibrium position in the "z" direction. 
It was first confirmed that this "rattle" was solely in the 
"z" direction (i.e., no significant "x" or "y" displacement) 
by extending the computation from 30 to 200 timesteps. It 
was shown that in a (100) direction in a body-centered cubic 
a narrow potential energy minimum was observed at (100) 
planes. The region of the (110) line between the octahedral 
void and the reference lattice site was contained within 
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this potential energy minimum. With this information as 
justification, it was determined that the velocity of the 
interstitial in the " z" direction could be completely damped 
at the end of each timestep to hasten the determination of 
the true equilibrium position. This technique was useful 
in restricting computer run time whenever an obvious "rattle" 
in the "z" direction developed. 

2 . The Dynamic Program 

The initial step of the dynamic program was the re- 
creation in the computer of the tungsten crystal, including 
the interstitial, after relaxation of the crystal had taken 
place and the interstitial had come to an equilibrium posi- 
tion. This was accomplished by utilizing the output of the 
static program as the input to the dynamic program. 

Before continuing with the problem, it was necessary 
to realize that the minimum energy required for the inter- 
stitial to escape from the crystal would be a function of 
the path of escape. To account for this "direction depend- 
ence", impact points were chosen in the surface plane of the 
crystal perpendicular to the direction of escape. The nega- 
tive "y" direction was always used as the direction of 
escape. Energy was then imparted to the interstitial which 
was, in turn, "aimed" at a specific impact point. The inter- 
stitial was then allowed to travel through the crystal using 
a timestep procedure based on a constant DTI. (The timestep 
procedure in the dynamic program is simpler than in the 
static program since energies are dominant throughout.) Each 
impact point (i.e., direction of escape) was subsequently 
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tested in the same manner using the same initial energy. 

If the interstitial escaped from the crystal in any di- 
rection tested, it was assumed that the initial energy 
applied to the interstitial was greater than the binding 
energy for that particular ion in that particular equili- 
brium position. In this case, the initial energy was de- 
creased and another survey of the impact points was taken. 

When the minimum initial energy which still allowed the 
interstitial to escape was determined, the binding energy 
for that particular ion in that particular equilibrium 
position was known. 

3 . Equilibrium Sites Viewed As Potential Wells 

One means of modeling equilibrium positions of 
interstitials in a crystal lattice is to consider each pos- 
sible equilibrium position as a three dimensional potential 
well. A foreign atom migrating through the lattice which 
reaches the confines of one of these potential wells with 
sufficiently small energy would fall into the potential well 
and remain there. In conjunction with this research, the 
character of these potential wells was also investigated. 

Two different approaches were taken to investigate 
the potential well aspect of the problem. First, the static 
program was used to investigate behavior of interstitials 
implanted at various positions around a previously deter- 
mined equilibrium position in a perfect lattice to determine 
which of these positions sought equilibrium. The second 
method offset the interstitial from its equilibrium position 
in the relaxed crystal and then allowed the crystal to undergo 



29 



further relaxation to determine whether the offset inter- 
stitial would seek the original equilibrium position. 

4 . Determination Of Possible Interstitial Sites 

In a body-centered cubic material there are twelve 
possible locations for a (110) split interstitial. In an 
infinite crystal these sites are equivalent; that is, no one 
site can be distinguished from another. When a finite crystal 
is considered, the presence of a crystal surface allows 
identification of three distinct interstitial sites. An "A" 
site is located in a (110) plane which contains the lattice 
site about which the split occurs and is closer to the 
crystal surface than the shared lattice site. A "C" site 
is also in a (110) plane containing the shared lattice site, 
but is located below the shared lattice site. A "B" site 
is located in the same (100) plane parallel to the surface 
that contains the shared lattice site, 
a. Implantation Procedure 

In the static program the interstitial was 
initially implanted in an offset position in the direction 
of the expected equilibrium position. This procedure al- 
lowed a smooth movement toward the suspected equilibrium 
position in a minimum number of timesteps. For the heavier 
interstitial atoms (xenon and tungsten) the lattice atom 
was also offset in its direction of suspected movement as 
a result of the presence of the large interstitial. 



30 



IV. PRESENTATION OF DATA 



A. STATIC SIMULATION 

1 . Initial Simulations 

After the minimum DTI procedure had been incorporated 
into the DTI decrementing process, investigations of tungsten 
interstitials in a tungsten crystal were made. In particular, 
a comparison between tungsten movement under the influence 
of an attractive potential and tungsten movement under the 
influence of a repulsive potential was sought. Equilibrium 
positions determined separately with these two potentials 
gave agreement within 0.02 lattice unit. 

Investigations were then extended to verify equili- 
brium positions near lattice sites 89 and 64 which had pre- 
viously been obtained by Tankovich [3]. Positions for argon, 
neon, krypton, and xenon were examined for each site, and 
results agreed with Tankovich' s data within 0.02 lattice 
unit, see Table I. As the mass of the interstitial atom 
decreased, its equilibrium position moved from the shared 
lattice site location previously reported [7,8] toward the 
center of the octahedral void. 

2. Site 39 

Investigations were then directed toward determining 
equilibrium positions of a split interstitial near lattice 
site 39, which is one layer below the surface of the crystal. 
Calculations were made separately for argon, neon, krypton, 
and xenon interstitials. At the end of thirty timesteps. 
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TABLE I 



Interstitial and Lattice Atom Displacements from Reference Site 



Interstitial 


Interstitial 
• Site 


Interstitial 

Displacement 


Lattice Atom 
Displacement 


Ax 


Ay 


Ax 


Ay 


Neon 


89C 


-0.86 


+0.87 


+0.06 


-0.06 


Argon 


89C 


-0.80 


+0.80 


+0.19 


-0.19 


Krypton 


89C 


-0.77 


+0.77 


+0.24 


-0.24 


Xenon 


89C 


-0.48 


+0.48 


+0.41 


-0.42 


Tungsten 


89C 


-0.38 


+0.39 


+0.35 


-0.35 
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all atoms of the crystal in each case had kinetic energies 
below thermal, which indicated that an equilibrium had been 
reached. An observation was made, however, that tended to 
abrogate this determination. The lighter interstitials, 
neon and argon, exhibited a definite affinity for the sur- 
face of the crystal. In all three possible sites (A,B, and 
C) , argon and neon interstitials were found to seek posi- 
tions significantly (~ 0.1 lattice unit) closer to the sur- 
face than the expected (110) split. Final positions for 
argon and neon interstitials in site 39A were located less 
than three tenths of a lattice unit below the surface of 
the crystal. Similar behavior, but to a lesser degree, 
was also observed with the heavier interstitials, krypton 
and xenon. Additionally, although the kinetic energy cri- 
teria indicated that an equilibrium had been reached, the 
small velocities that were available at the end of each 
timestep were in such a direction as to allow escape from 
the crystal if these velocities could be maintained. In 
fact, although velocities were halved at the end of each 
timestep, a comparison of velocities over the last five 
timesteps of the computer run showed that the negative "y" 
(direction of escape) velocity component of the interstitial 
at the end of timestep thirty was actually only twenty per- 
cent lower than the same velocity component at the end of 
timestep twenty five. In other words, even with fifty per- 
cent damping applied during each timestep, velocities actu- 
ally decreased by only twenty percent over 5 _ timesteps. It 
was also observed that during these last five timesteps the 
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maximum movement allowed by the most active atom of the 
crystal (the DTI) had reached its minimum value, 0.0005 
lattice unit. These observations suggested that the small 
movement allowed during these timesteps combined with the 
damping factor might be "forcing" the interstitial to ex- 
hibit equilibrium criteria. It was concluded that equili- 
brium positions near site 39 were unstable, if they exist 
at all. 

3. Site 14C 

Until these investigations little thought was given 
to the possibility of an equilibrium position in site 14C. 
With no precise knowledge available concerning the actual 
quantitative value of the damping experienced by a foreign 
interstitial implanted in a crystal lattice and with the 
possibility of equilibrium sites near lattice site 39, 
some credence had to be given to the possibility of equili- 
brium positions near surface atoms of the crystal. It was, 
therefore, decided to investigate the possibility of an 
equilibrium position in site 14C. It was postulated that 
in the case of light atoms (neon, for example) the inter- 
stitial position for the 14C site should be deeper in the 
crystal than the 39A site, as a result of the greater re- 
pulsion of the lattice atom with which the site was shared. 
Simulations showed that the 14C interstitial site was located 
at a distance of 0.43 lattice unit below the surface of the 
crystal while the 39A interstitial site was located at a 
distance of 0.26 lattice unit below the crystal surface. 

When the displacement of this interstitial site was compared 
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with the 64C interstitial site (0.84 lattice unit below 
the next higher lattice plane) , the effect of the surface 
of the crystal on interstitials in close proximity to the 
surface was clearly exhibited. Velocity characteristics 
similar to those of site 39A were also observed indicating 
that the amount of damping applied and the DTI could be 
forcing the interstitial to exhibit equilibrium properties 
in this site. In general, it can be said that interstitial 
equilibrium sites in the first two layers of the crystal 
are ill defined in position if they exist at all. 

4 . Force Calculations 

As mentioned previously, in the course of these 
investigations it was determined that energies nearly al- 
ways dominated forces for the timestep ranges used in the 
calculations, and, consequently, the timestep was nearly 
always chosen as a function of the energy of the most ener- 
getic atom. To ascertain whether the use of timesteps deter- 
mined by maximum forces would yield better or significantly 
different results, the static program was modified so that 
the timesteps were determined strictly as a function of the 
maximum force. This was accomplished by first printing out 
the maximum force observed during each timestep. By sur- 
veying these maximum forces, several values of force were 
chosen to be used as test forces. As the program was sub- 
sequently run, the maximum force in each timestep was compared 
with these test forces, and a DTI for that timestep was as- 
signed based on the results of that comparison. This 
assigned value of DTI was then used to compute the next 



35 



timestep interval. These calculations of equilibrium 
positions based on forces agreed with previous energy cal- 
culations of equilibrium positions within 0.03 lattice 
unit. 

It should be realized here that the terminology 
"force calculations" and "energy calculations" do not imply 
any significant difference in method for determining equili- 
brium positions. They are merely two different procedures 
for determining the maximum displacement which will be al- 
lowed by - the most energetic atom of the crystal during the 
next timestep. The importance of this parameter (DTI) is 
in the effect it has in insuring the smooth movement of the 
interstitial to its equilibrium position. 

5 . Investigation of Potential Wells 

Viewing interstitial equilibrium positions as po- 
tential wells of varying depths is a convenient means of 
modeling the entrapment of foreign atoms by lattice struc- 
tures. In order to minimize the effects of the crystal 
surface on investigations, interstitial site 89C, which had 
exhibited good stability and almost perfect (110) splitting 
in previous testing, was chosen for investigation of the 
characteristics of interstitial potential wells. 

a. Potential Well Studies in a Perfect Lattice 

The first approach to the study of the inter- 
stitial potential well utilized essentially the same method 
that had been originally employed to locate equilibrium posi- 
tions. It was postulated that any interstitial atom that was 
implanted in a perfect lattice at or near the coordinates of 
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the equilibrium position which had been previously deter- 
mined for interstitial site 89C would seek the same equili- 
brium site. Interstitial atoms could be implanted further 
and further from the previously determined equilibrium 
position until they no longer returned to it. In this man- 
ner a map of the size of the potential well around site 89C 
could be obtained. 

The results obtained in investigations of xenon 
in site 89C are presented in Figure 3. In this figure each 
coordinate intersection represents an implantation site 
which was tested. The arrows drawn from these implantation 
sites indicate the direction of movement of the interstitial 
from that specific implantation site. The tip of the arrow 
represents the interstitial position after thirty timesteps. 

In considering the data obtained in these in- 
vestigations, two significant observations were made. First, 
an interstitial implanted at the previously determined equili 
brium position (coordinates (4.52, 3.48) in Figure 3 ) moved 
from this position to another position which also met the 
criteria for equilibrium. Secondly, all other implantation 
sites also moved to positions which met equilibrium criteria. 
An analysis of the program print out data provided informa- 
tion of secondary importance. All implanted atoms exhibited 
an initial movement in the general (110) direction even 
though final positions were not necessarily in that direction 
Similar investigations were conducted using argon and neon 
interstitials with similar results. 
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In light of this unexpected behavior in the 
perfect lattice the stability of a foreign atom, inserted 
into the crystal as a replacement for the lattice atom in 
site 89, was examined. Investigations were made using neon, 
argon, krypton, and xenon as the replacement atoms. The 
data from these runs are tabulated below. 



TABLE II 

Results of Replacement Runs in Site 89 



Replacement 


Y Displacement (L.U.) 


Final Kinetic 


Final Potential 


Atom 


After 30 Timesteps 


Energy (eV) 


Energy (eV) 


Neon 


-0.049 


0.0000 


1.0444 


Argon 


-0.036 


0.0008 


3.5118 


Krypton 


-0.029 


0.0001 


3.5120 


Xenon 


-0.020 


0.0000 


1.3297 



The y displacement values are sufficiently small 
that site 89 must be presumed to be a stable replacement site. 
This conclusion is further confirmed by observing that the 
potential energies are also significantly lower than those 
obtained for the same atoms when acting as interstitials, 
see Table III. 



TABLE III 



Interstitial Potential 

Interstitial 

Atom 

Neon 

Argon 

Krypton 

Xenon 



Energies for Site 89C 

Final Potential 
Energy (eV) 

4.7 

16.3 

15.5 

9.1 
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To further check the stability of this site, 
a kryptron replacement atom was then initially offset from 
the site 89 coordinates, 0.7 lattice unit in the "x" and 
"y" directions, and the program was run again. The krypton 
atom moved precisely along the (110) line back toward site 
89. At the end of thiry timesteps, the krypton atom was 
located 0.4 lattice unit in the "x" and "y" directions, 

"x" and "y" velocity components were still in the direction 
of site 89, and potential energy had decreased from 54 eV 
in the iiiitial offset position to 8.7 eV at the end of 
thirty timesteps. 

(1 ) Interpretation of Results 

Interpretation of the results reported in 
the previous section led to a re-examination of the concept 
of the equilibrium positions of interstitials and raised a 
question concerning the validity of solely using relaxed 
crystals for determinations of equilibria. 

The equilibrium positions obtained by using 
a perfect crystal were a function of the implantation site. 
This dependency of the final equilibrium position on the 
implantation site suggested that the movement of the inter- 
stitial to an equilibrium site can not be explained in terms 
of the lattice structure "forcing" the interstitial to its 
lowest energy configuration. Rather, these results sug- 
gested that the actual mechanism is a process of local 
"liquefaction" of the lattice as it adjusts to the presence 
of the interstitial combined with interstitial movement caused 
by interaction with the lattice atoms. This, in turn. 
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suggested that a more accurate means of localizing equili- 
brium positions might be to utilize a crystal which had 
already been allowed to "relax", or adjust to the presence 
of the interstitial. (See next section.) 

The movement of atoms from nine different 
implantation sites in an area two tenths of a lattice unit 
per side to nine different positions which exhibit equili- 
brium criteria indicated that equilibrium positions are not 
precise positions coordinate - wise within the time range 
of these computations. 

b. Potential Well Studies in a Relaxed Lattice 

The conclusions reached as a result of studies 
of potential wells in a perfect crystal prompted similar 
studies in a crystal which had been allowed to adjust to 
the presence of the interstitial. An argon interstitial 
and interstitial site 89C were chosen for investigation. 

The relaxed lattice was obtained by using the results of the 
static simulation of argon in site 89C which had first in- 
dicated that an equilibrium position had been reached. These 
results were read into the computer as initial positions for 
the lattice atoms and the interstitial. The interstitial 
was then offset from its position and a new static simu- 
lation was performed. An array of sites was chosen for in- 
vestigation in this manner. The results of these simulations 
with various interstitial offsets are shown in Figure 4. 

The representation in Figure 4 is analogous to the repre- 
sentation in Figure 3. The start point for each run is 
numbered (1-9 and A,B). 
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The interstitial movement depicted in Figure 4 
was encouraging in many respects. An argon interstitial 
placed successively at each initial site (numbered 1-9 in 
Figure 4) appeared to head for roughly the same region of 
the crystal. The potential energy of the interstitial at 
the conclusion of each run ranged in value from 8.8 to 9.8 
eV. A comparison of these energies with the energy of the 
interstitial in the initial equilibrium position found in 
the perfect lattice (~16 eV) indicated that this equilibrium 
region is much more stable than the position determined pre- 
viously. While not precisely defined in position, this 
equilibrium region was located in roughly the (110) di- 
rection from the shared lattice site. 

Since this equilibrium region appeared to be 
located nearer to the shared lattice site (site 89) than 
offset start points 1-9, two additional computer runs (labeled 
A and B in Figure 4) were made with start points bordering 
this equilibrium region. These runs indicated that, even 
in a relaxed crystal, the equilibrium reached was still 
somewhat a function of the initial position of the inter- 
stitial. In both of these runs, equilibrium was reached by 
movement in the (110) direction, and the potential energy of 
the interstitial after thirty timesteps (~ 10 eV) was lower 
than the potential energy originally determined in the per- 
fect lattice. 

(1) Interpretation of Results 

The results obtained in the relaxed lattice 
seemed to reinforce the conclusions drawn during the study 
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of potential wells in the perfect lattice. Basically, the 
relaxed crystal studies gave a picture of equilibrium posi- 
tions much more in tune with what might reasonably be ex- 
pected in nature. It seems clear that the equilibrium site 
of an interstitial is an ill-defined region in the general 
(110) direction from the shared lattice site. Determination 
of this equilibrium region cannot be made based solely on 
the kinetic energy of the interstitial; potential energies 
must also be considered. 

. This equilibrium region could have b een 

postulated from a consideration of the rms displacement of 
a tungsten atom in a tungsten lattice. Houska [18] has 
measured this rms displacement to be 0.049 lattice unit at 
300°K. The ordinary thermal activity of the lattice atoms, 
then, can be expected to cause the equilibrium positions of 
interstitials to fluctuate over some region. This equili- 
brium region might best be described as roughly a cylinder 
whose axis is in the <110) direction and whose height and 
radius are some function of the relationship between inter- 
stitial mass, lattice atom mass, and the interatomic poten- 
tial function. The equilibrium positions determined in these 
investigations appeared to be centered in the <110) direction 
at the intersection of the (110) and (100) planes through 
the lattice atom and covered a section of the <110) line 
on the order of 0.2 of a lattice unit long. This seems to 
be a reasonable range of equilibrium positions for the re- 
latively light argon interstitial in a tungsten crystal. 
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The results of these investigations have 
indicated the usefulness of conducting implantation studies 
in relaxed crystals. The movement of interstitials to equi- 
librium in the relaxed crystal was much less dependent upon 
initial positions, and all sites investigated showed a pre- 
ference for positions in the <110) direction. Such was not 
the case in perfect lattice studies. 

B. DYNAMIC SIMULATIONS 

Concurrent with the investigations into the potential 
wells of the equilibrium sites, dynamic simulations were 
conducted using the argon interstitial in site 64A. Since 
Tankovich [3] had investigated an argon interstitial in 
site 89C and had determined that the binding energy of this 
site was in the range 2-4 eV, 4 eV was chosen as the initial 
energy for the dynamic tests. It was determined that the 
binding energy of the argon interstitial in site 64A was 
between 0.02 and 0.04 eV. At 0.04 eV escape of the inter- 
stitial was noted for several of the twelve impact points 
tested. At 0.02 eV the interstitial first moved slightly 
toward the surface- of the crystal, then changed directions 
and moved to a position deep inside the crystal while gaining 
considerable energy. This phenomena had been seen and in- 
terpreted previously as the expected behavior of a stable 
site in which the interstitial has been required to move 
too large a distance in one timestep. Since the timestep is 
constant (0.1 lattice unit) in the dynamic simulation, an 
interstitial that is oscillating in a stable configuration 
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can be required to move out of its stable position. When 
this movement is directed into the crystal, it can be assumed 
that the original position was stable . 

This binding energy falls in the range of the first 
desorption peak measured by Kornelsen and Sinha (see Figure 
5 in Ref. 1) . This result indicated that interstitial sites 
formed with lattice atoms in the third crystal plane (cor- 
responding to the y = 2 plane of these calculations) , such 
that the interstitial site was in the (110) direction above 
the lattice atom, are the sites nearest the surface of the 
crystal stable enough to entrap measurable amounts of argon. 
This excellent concurrence with experiment further sub- 
stantiated the conclusion drawn earlier in these investi- 
gations about the doubtful character of positions near sites 
39 and 14 which had exhibited equilibrium criteria. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



During the course of these investigations the possibility 
that the DTI decrementing procedure and the procedure for 
incorporating damping into the problem might be forcing in- 
terstitials to exhibit equilibrium properties arose. At 
first, the DTI was thought to possibly be restricting atom 
movement to such an extent during the last few timesteps of 
a calculation (i.e., when DTI had reached its minimum value 
-0.0005 lattice unit) that atom movement and, consequently, 
kinetic energy could no longer be used as equilibrium cri- 
teria. This became more significant when the small DTI and 
the damping factor were considered together. The damping 
factor seemed particularly appropriate for scrutiny near the 
surface of the crystal and when DTI was small, since the pro- 
bability of collisions and close approaches and the subsequent 
damping of atom motion could be expected to decrease in both 
instances. The excellent agreement with experiment of the 
results of the dynamic and static simulations suggests, how- 
ever, that the DTI and damping procedures employed might be 
adequate for these simulations except when atom movement is 
particularly close to the surface of the crystal. It is 
recommended that future investigations explore the use of a 
damping factor which varies with decreasing DTI and movement 
toward the surface of the crystal. 
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These investigations have shown the value of using a 
relaxed crystal in the computer simulation of lattice dynamics. 
The' relaxed crystal simulation indicated that: 

(1) The equilibrium position initially sought by an 
interstitial is a function of implantation position. 

(2) The mechanism associated with the establishment of 
an equilibrium position seems to be a combination of inter- 
atomic interaction and local liquefaction of the crystal 
structure in the vicinity of the interstitial. 

(3) The character of equilibrium potential wells is more 
readily observable in a relaxed crystal. 

(4) The actual equilibrium position of an interstitial 
seems to be some region in the general (110) direction from 
the reference (shared) lattice site. 

Again, however, for determining binding energies, the 
original procedure for determining equilibrium positions may 
yield adequate results. Such was the case for the dynamic 
simulations reported here. It is recommended that replace- 
ment interstitials and other split interstitial sites undergo 
dynamic simulation in an attempt to correlate other desorption 
peaks as determined by Kornelsen and Sinha [1] with inter- 
stitial sites. 
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APPENDIX A 



SUBROUTINE PLUCK 

Subroutine Pluck was developed to allow the use of a 
crystal smaller than the original model in dynamic testing. 
This subroutine was developed so that the only pieces of 
information required as input to the subroutine were the 
number of crystal planes desired in the "x" and "z" 
directions and the number of the lattice site desired at 
the center of the PLUCK crystal. Once the size of the 
PLUCK crystal is decided upon, the subroutine stores the 
original numbers of the atoms in the PLUCK crystal in an 
array. This allows reference to any atom by its original 
number throughout computation. The atoms of the PLUCK 
crystal were numbered consecutively, and the number of atoms 
in the PLUCK crystal was assigned a variable name. In this 
manner, a minimum of adjustment was required in the dynamic 
program when the PLUCK crystal was used. 

SUBROUTINE PLUCK is included in this appendix in its 
most general form. Parameters are listed below for two 
different sizes of PLUCK crystals. Sizes refer to the 
number of "x" and "z" planes in the PLUCK crystal. "Y" planes 
from the surface of the crystal through the two planes below 
the plane of the lattice site under investigation are always 
included. 
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PLUCK PARAMETERS FOR DIFFERENT SIZED CRYSTALS 



(crystal sizes refer to number of planes of atoms in x and 
z directions) 



Parameter 

IXNEW 

IZNEW 

NM1 

Nil 

NX1 

NII31 

NII41 

NIINCl 

IIINCl 

NMINCl 

NXINC1 

NMINC2 

NX 2 

NM2 

NI2 

NII32 

NII42 

NIINC3 

I I INC 2 

NMINC3 

NXINC2 

NMINC4 



7x7 Crystal 
7 

7 
5 

8 

16 

10 

5 

1 

4 

4 

9 

-1 

9 

4 

8 

4 

11 

2 

3 

3 

16 

1 



5x5 Crystal 
5 
5 

3 

14 

4 

16 

9 

3 
2 
2 
9 
1 
9 

4 
8 

10 

15 
2 
3 

3 

4 

-1 
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•SUBROUTINE PLUCK' 

aftc fc* ## A#######:*;## 



COMMON/ COMl/RX (500) ,RY(500) , RZ(500) , LCUT(500), 
1LL,LD, I TYPE,NVAC 

C0MM0N/CCM4/ IX, I Y, I Z , SCX , SCY , S C Z , IDEEP,D1X,D1Y,D1Z, 
1 IVACX.I VACY,IVACZ 
COMMON /COM 1 0/ I XNEW , IYNEW, I ZNEW, I I 

COMMON/ C0M11/RXNEW I (250) ,RYNEWI ( 250 ) , RZNEW I ( 250) , 
1KEEP ( 250) , NNUM ( 2 50 ) 

1500 IXNEW=IXNEW 

IYNEW=I VACY+3 
IYNEW=I YNEW 
NM=NM1 
NI =N 1 1 
11 = 2 
MM=0 
NX=NX1 
N I 1 3 = N I I 31 
NI 14 =N I 141 

I F ( I YNEW.EO. 3) GO TO 1514 
I Ft IYNEW .EG. 5) GO TO 1514 
1505 DO 1509 1=1 I,NM 
NNUM ( I ) =NI 
1509 N I = N I + 1 

N I =N I +N I I NCI 
I I = I I + I IINC1 
NM=NM+NM INC 1 

IF { I I . LE. NX) GO TO 1505 
N X=N X+ N X I NC 1 
NI=NI+NI 14 
MM=MM+ 1 

I F ( I YNEW.EO. MM ) GO TO 1600 
NM=NM+NM INC2 
GO TO 1515 

1514 NX=N X2 
NM=NM2 
N I =NI 2 

NI I 3 = N I 132 
NI I4=N I 142 

1515 DO 1520 I = I I , NM 
NNUM ( I ) = N I 

1520 N I = N I + 1 

NI =NI+NI INC3 

I 1= I I + I I I NC 2 

NM=NM+NM I NC3 

I F ( I I . LE. NX) 30 TO 1515 

NX=NX + NX INC 2 

N I =N I +NI 13 

MM=MM+ 1 

I F ( IYNEW. EO. MM ) GO TO 1600 
NM = NM+ NMI NC4 
GO TO 1505 
1600 11=11-1 

RXNEWI ( 1 ) = RX ( 1 ) 

RYNEWI ( 1 ) =RY ( 1 ) 

RZNEW I ( 1 ) = R Z ( 1 ) 

KEEP ( 1 ) =1 
NNUM ( 1 ) = 1 

1700 DO 1750 1=2,11 

RXNEWI ( I ) =RX( NNUM ( I ) ) 

RYNEWI ( I ) = RY ( NNUM( I ) ) 

RZNEW I ( I } = RZ( NNUM( I ) ) 

KEEP ( I ) =NNUM( I ) 

1750 CONTINUE 
RETURN 
END 
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APPENDIX B 





COMPUTER PROGRAM GLOSSARY 


ALPHA: 


Input Morse potential parameter. 


BSAVE: 


Target Mass/ (target mass + bullet mass) ; 
distributes potential energy between target and 
bullet. 


BIND: 


Negative of the total potential energy (TPOT) 
at time zero. 


BMAS : 


Mass of bullet in AMU. 


BULLET : 


Alpha-numeric array for point defect material. 



CFO, CFl , CF2: Force parameters of cubic fit between Morse 





and Born-Mayer functions. 


CGB1, CGB2 : 


: Morse potential parameters. 


CGDl , CGD2 : 


: Morse potential parameters. 


CGF1 , CGF2 : 


: Morse force parameters. 



CPO, CPI, CP2, CP3: Potential parameters of cubic fit 

between Morse and Born-Mayer functions. 

CVD: CVR X 10 Converts lattice units to meters. 

-19 

CVE: 1.6 X 10 : Converts electron volts to Joules. 

CVED: CVE/CVD: A ratio used to avoid repeated division. 

-27 

CVM: 1.672 X 10 : Converts atomic mass units to kilograms. 

CVR: LU in angstroms; converts lattice units to angstrom 
units . 

DIX, DIY, DIZ: Displacement coordinates for location of 



DCON : 


interstitial from reference atom, NVAC. 
Input Morse potential parameter. 
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DDTIF : 


The minimum value that DTI is allowed to assume. 


DDTII : 


The initial decrement of DTI. 


DFF : 


ROE-DIST. The distance closer than ROE that an 
atom is to the primary. 


DIST: 


Distance between any two atoms. 



DLPE: TLPE-TLPE0: The change in total local potential energy 
since time zero. 

DRX, DRY, DRZ : x,y,z components of DIST. 



DT: 


Length of a timestep in seconds. 


DTE, DTEl : 


The two possible alternatives of the timestep 
computed from maximum energies. 


DTF, DTF1 : 


The two possible alternatives of the timestep 
computed from maximum forces. 


DTI: 


Number of lattice units most energetic atom may 
move in one timestep. 



DTI (I), DT2 (I), DT3 (I), DT4 (I): Vector arrays which save 



DTSTEP (I) 


the possible choices of timestep determined in 
the "energy" method. 

: Vector array which saves the timestep interval 
chosen for each timestep. 



DTOD: DT/CVD -- A ratio used to avoid repeated division. 
DTOM :DT/PTMAS — A ratio used to avoid repeated division. 
DX(I) , DY(I), DZ(I): Change in position of ith atom from 





initial position at time zero. 


EMAX: 


The maximum energy encountered in any cycle. 


EV: 


Primary energy in electron volts. 


EVR: 


Primary energy in kiloelectron volts. 
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EXA, EXB : 



EXA, E> 


□3: Input Born-Mayer potential function parameters 

for the target. 


F2: • 


Square of the force on a specific atom. 


FA: 


The component force increment on an atom. 



FDTI : DTI X CVD: A parameter used to determine DT by 





maximum energy method. 


FM: 


A small number used in checking potential energy 
zero point. 


FM2: 


FM squared. 


FMAX: 


Maximum total force on the most stressed atom 
in the crystal. 



FOD: FORCE/DIST: A ratio used to avoid repeated division. 



FORCE : 


Numerical value of the force function with a 
variable parameter. 


FORMAX 


(I) : Vector array which saves the maximum force in 
each timestep. 



FX(I), FY(I), FZ(I): x,y,z, components of total force on 



FAX: 


an atom. 

Born-Mayer force function parameter. 



HBMAS: % BMAS : A ratio used to avoid repeated division. 

HDTOD: ^ DTOD: A ratio used to avoid repeated division. 

HDTOM: \ DTOM: A ratio used to avoid repeated division. 

HDTOMB: % DTOMB : A ratio used to avoid repeated division. 
HTMAS: ^ TMAS : A ratio used to avoid repeated division. 



11: 


Variable in cubic fit subroutine. 


13: 


Variable in cubic fit subroutine. 


IDEEP: 


First fixed layer. 
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IFMXAM 


(I) : Vector array which saves the atom number which 
experiences the maximum force. 


IH1 : 


Alpha numeric array for program title. 


IH2 : 


Alpha numeric array for Morse function parameters 


IHB : 


Alpha numeric array for bullet element. 


IHS : 


Alpha numeric array for type and orientation of 
crystal. 


IHT: 


Alpha numeric array for target element. 


II: 


Number of atoms in a crystal using subroutine 
PLUCK. 


I LAY: 


Number of free (mobile) layers. 


IN: 


Odd-even integer used to determine atom site 
establishment. 


IP: 


Subscript value of atom. Void in subroutines 
STEP and ENERGY. 


IQ: 


Parameter that determines whether or not a self 
defect is to be given a repulsive potential or 
a composite attractive - repulsive potential. 


ISHUT: 


A parameter used to shut down the program. 


IT: 


Unsealed fixed point x coordinate used in 
lattice generation. 


ITT: 


Odd-even integer used to determine atom site 
establishment. 


I TYPE: 


Parameter used to determine the type of point 
defect: vacancy, interstitial, or replacement. 



IX, IY, IZ: Number of x,y,z planes of cyrstal. 

IXNEW, IYNEW, IZNEW: Number of x,y, and z planes in the 
PLUCK crystal. 
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J2: 


Variable in the cubic fit subroutine. 


JJ: 


Parameter in the BCC(lll) lattice generation 
subroutine. 


JT: 


Unsealed y coordinate used in crystal generation. 


JTS : 


Variable used to establish atom sites. 


JTT: 


Variable used to establish atom sites. 


KEEP (I) : 


Vector array which saves the original atom 
numbers of the PLUCK crystal. 


KF: 


Final K in LOCAT (K) assigned to an atom. 


KT: 


Unsealed z coordinate used to establish atom 
site . 


LCUT (I) : 


Used to identify an ith atom which is not 
included in calculations. 


LD: 


The highest numbered atom in the mobile layers. 


LL: 


The highest numbered atom in the entire crystal. 


LOCAT (K) : 


Dimensioned variable that remembers the numbers 
of the atoms within a radius ROEL of the primary 
at time zero. 


LS: 


Variable associated with each of the nine lattice 
generator subroutines. 


MCRO: 


One number higher than the order of the fit 
between the Born-Mayer and Morse potentials, 
always 4 in this simulation. 


ND: 


Data output increment, in numbers of timesteps. 


NEW: 


Parameter used to determine whether or not atom 
numbers have been stored in LOCAT (K) . 


NNUM (I) : 


Vector array used in Subroutine PLUCK to re- 
number atoms . 
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NPAGE : 



NRUN: 



NS: 

NT: 

NTT: 

NVAC: 



PAC: 

PBMAS : 

PEXA, PEXB : 

PFPTC : 

PFXA: 

PKE (I) : 
PLANE: 

POT: 

PPE (I): 
PPEINT (I) : 

PPENCK: 

PPENEG (I) : 

PPEPCK: 

PPEPOS (I) : 



Page numbering variable. 

Parameter used to determine whether or not to 
read additional data cards. 

Initial print statement timestep number. 

Timestep number. 

Timestep number limit before shutdown. 

An atom number used to establish point defects 
or used as a reference point for interstitial 
placement. 

Parameter for bullet force function correction. 
Primary mass in kilograms. 

Input Born-Mayer potential function parameters 
for the bullet-target interaction. 

Primary force function evaluated at ROE. 

Primary force function parameter. 

Kinetic energy of the ith atom. 

Alpha-numeric array for lattice orientation. 
Potential energy between two atoms. 

Potential energy of the ith atom. 

Vector array that saves the difference in 
potential energy before and after implantation. 
Potential energy check value which determines 
potential energy decreases which will be printed 
Vector array which saves potential energy 
decreases. 

Potential energy check value which determines 

which potential energy increases will be printed 

Vector array which saves potential energy in- 
creases. 
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PPESAV (I) : Vector array which saves the initial potential 



energy of lattice atoms before implantation. 
PPKEEP (I) : Vector array which saves potential energy dif- 



PPTC : 


ferences between perfect crystal and relaxed 
crystal . 

Primary potential function evaluated at ROE. 



PTE (I) : Total energy of the ith atom (potential + 





kinetic) . 


PTMAS: 


Target mass in kilograms. 


RE: 


Input Morse potential parameter. 


RO: 


Spacing constant in BCC(llO) lattice generation 
subroutine. 


ROE: 


Nearest neighbor distance. 


ROE 2 ; 


ROE squared. 


ROEA: 


Maximum cut-off for Born-Mayer potential. 


ROEB : 


Minimum cut-off for Morse potential. 


ROEC :’ 


Maximum cut-off for Morse potential. 


R0EC2: 


ROEC squared. 


ROEL: 


Radius inside of which local potential energy 
is found. 


R0EL2 : 


ROEL squared. 


ROEM: 


ROE— DTI , Region in which modification of 
repulsive force must be made. 



RX(I), RY(I), RZ(I): x,y,z, coordinates of an ith atom 
at any time. 

RXI (I) , RYK(I), RXI(I): x,y,z, coordinates of an ith atom's 
initial position. 
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RXK(I) , 


RYK(I), RZK(I): x,y,z coordinates of temporary 
position of an ith atom durinq force cycle. 


RXNEWI, 


RYNEWI, RZNEWI : Vector arrays which contain the 
x,y,z, coordinates of the atoms of the PLUCK 
crystal. 


RXSAVE , 


RYSAVE, RZSAVE: x,y,z, coordinates of the impact 
point in the dynamic progran . 


SAVE: 


h POT. 



SCX, SCY, SCZ: x,y,z, coordinate scale factors. 



SSCZ : 


• A z scale factor used for the FCC (111) lattice 

generator subroutine. 


START: 


An optional timing variable, not used in this 
simulation. 


SUM: 


Variable in cubic fit subroutine. 


TARGET : 


Alpha-numeric array for target material. 


TSAVE : 


Bullet mass/ (target mass+bullet mass) ; 
distributes potential energy between target 
and bullet. 


TE: 


Total energy of all crystal atoms (kinetic + 
potential) . 


TEMP: 


Temperature of lattice in degrees K elvin not 
used in this simulation. 


TFAC : 


A time factor ratio used to determine DT by 
maximum force method. 


TFACB : 


TFAC for the bullet. 


THERM: 


Thermal energy of atom not used in this simu- 
lation . 


TIME: 


Elapsed problem time in seconds. 
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TLPE : 


Total local potential energy of atoms within 
a radius ROEL. 


TLPE0 : 


TLPE at time zero. 


TMAS : 


Target atom mass in AMU. 


TPKE : 


Total kinetic energy of all crystal atoms. 


TPOT: 


Total potential energy of all crystal atoms. 


VSS : 


Storage variable for velocity components. 


VS(1) , 


VY(1), VZ(1): x,y,z components of ith atoms velocity 



X, Y, Z: Unsealed coordinates used in crystal generation 

XNVAC, YNVAC, ZNVAC : The initial displacement (in LU) of 
atom NVAC. 

YLAX(I): Relaxation in -y direction of ith layer in L.U. 

ZP: Floating point form of JTT. 
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APPENDIX C 



ERROR IN THE LITERATURE 

It was discovered during the course of these investi- 
gations that Equation (20) of Ref. 10 (corresponding to 
equation 6A of the report) was incorrect. 

Equation (19) of Ref. 10 is, 

AT A AX ± /v i 

where AT is the timestep interval, A*^ is the displacement 
of the i th particle, and v i is the velocity of the i th 
particle . 

To express the timestep in terms of the energy of the 
particle with the maximum energy, Ax^ is defined as the 
displacement of the particle with the greatest energy and 
is replaced by the symbol D. T^ is this particle's energy, 
which is the largest kinetic energy at end of each timestep. 
If T is expressed as 



m 



T = % mv , 



m 



m' 



substitution into equation (19) yields, 

h 



AT = D/(2 T m /m) 2 , 



or 



vice 



AT = D (m/2 T m ) 






AT = D(2m/T m ) 



2 . (Equation (20) Ref. 10) 

It should be noted that this error has no affect on the 
calculations reported here. D (DTI in these calculations) 
is a parameter which has been specifically chosen in order 
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to provide a timestep interval over which there are no 
significant force changes. The factor of 2 which is intro- 
duced in this correction would merely have required a sub- 
sequent alteration of the DTI 1 s chosen so that the timestep 
interval would remain essentially the same. 
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Figure 1. Atom Numbering Procedure. 
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THE FUNCTION OF SUBROUTINE PLUCK 




Figure 2. The Function of Subroutine Pluck. 
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Figure 3. Xenon Behavior in a Perfect Lattice. 
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Figure 4. Argon Behavior in a Relaxed Lattice. 
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COMPUTER PROGRAM 
FOR STATIC SIMULATIONS 



STATIC SIMULATIONS WERE USED TO INVESTIGATE THE 
EQUILIBRIUM POSITIONS CF INTERSTITIAL ATOMS IMPLANTED IN 
A TUNGSTEN CRYSTAL. FOUR DIFFERENT CONFIGURATIONS OF THE 
STATIC PROGRAM WERE USED AS THE INVESTIGATIONS PROGRESSED. 

these Configurations were used to 

(1) DETERMINE EQUILIBRIUM POSITIONS BY 'ENERGY' 
CALCULATIONS. 

(2) DETERMINE EQUILIBRIUM POSITIONS BY 'FORCE' 
CALCULATI ONS, 

(3) PRINT OUT THE SMALLAR CRYSTAL DETERMINED BY 



SUBROUTINE PLUCK FOR USE IN INITIAL DYNAMIC 
TESTING, AND 

(4) INVESTIGATE THE POTENTIAL WELLS IN THE RELAXED 
CRYSTAL . 

THE 'ENERGY' CONFIGURATION OF THE STATIC PROGRAM IS PRESENT- 
ED BELOW WITH DIFFERENCES FROM THIS PROGRAM REQUIRED BY 
OTHER CONFIGURATIONS INCLUDED AND DISCUSSED AT APPROPRIATE 
POINTS. ADDITIONALLY, BRIEF COMMENTS ARE INCLUDED AT VARIOUS 
POINTS TO CALL ATTENTION TO SIGNIFICANT PROCESSES OF THE 
PROGRAM. 



C 

C DIMENSION VECTOR ARRAYS USED EXCLUSIVELY IN THE MAIN 

C PROGRAM. THE DIMENSIONING SCHEME FOR EACH CONFIG- 

C URATION IS GIVEN BELOW IN ITS ENTIRETY. 

C 






ifctk A* jfc # £ # ^ £ # # # * r$c * sft # # # 

C 'FORCE CONFIGURATION* 

# * # £ ## & # £ # Jfc ■>£ ^ 5$C & sk * A 3$C # 3fc XS 5$C # # 

DIMENSION VX(500),VY(500) ,VZ(500) , P KE ( 5 00 ) 

DIME NS I ON DX (500 ) , DY( 500) , DZ ( 500 ) ,PTE( 5 00) 

D I ME NS I ON PPESAV ( 500) , PPEINT ( 500 ) , PPEPOS (500 ) , 

1 PP ENEG ( 500 ) 

DIME NS I ON FSTACC ( 1 00 ) ,FSSACC(100 ) , FORMA X ( 100 ) , 

1 1 FMXAM ( 100) 

DIMENSION DTFSP(50 ) ,DTIASP( 50), DTESP( 50) ,DTIVSP( 50) , 
1DETI ( 50) 

DIMENSION P PKE E P ( 50 0) 

DIMENSION DTSTEP(200) 



C 






•OTHER CONFIGURATIONS’ 



DIMENSION VX(500) ,VY( 500) ,VZ( 500) ,PKE( 500) 
DIMENSION DX ( 5 00 ) , DY ( 5 00 ) , DZ ( 50 0 ) , P T E ( 5 0 0 ) 
DIMENSION PPESAV( 500) , PPEINT (500) , PPEPOS (500 ), 
1 PP ENEG ( 500 ) 

DIMENSION FSTACC (100) , FSSACCllOO ) ,FORMAX( 100), 
1 I FMXAM ( 100) 

DIMENSION DTI ( 50 ) , DT2 ( 5 0 ) , DT 3 ( 50 ) , DT4( 50) 
DIMENSION P PKE E P ( 5 00 ) 

DIMENSION DTSTEP(200) 



A* £ $ sjc## + * + 5{c # # 5je afe 
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PRESCRIBE COMMON STORAGE OF VARIABLES AND VARIABLE 
ARRAYS REQUIRED IN SUBROUTINES. 



C 



COMMON/ C0M1/RX ( 500 ) ,RY ( 500 ) , RZ( 500) , LCUT( 5 00) , 

ILL ,LD,ITYPE ,NVAC 

COMM ON /COM 2/ I H 1 ( 20) ,IH2(8),IHS(10),IHB(6) ,IHT(6) , 

I TARGET (4 ),TMAS , BULLET (4) , BMAS , PL AN E, TEMP , THERM, 

1DDTI I , DDTI F 

COMMON/COM3/RXI (500),RYI(500),RZI(500),CVR,EVR, 

1 NT, TI MF , DT , DTI , I LAY , RXK ( 500 ) , RYK ( 500 ) ,RZK{ 500) 
C0MM0N/C0M4/I X, I Y, IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z, 

I IVACX, IVACY, IV ACZ 

COMMON/CCM5/ROE, ROE2, ROEM, EX A, EXB, PEXA, PEXB, FXA, PFXA, 
1 IQ,T SAVE , BSAVE 

COMMON/ C0M6/FX ( 500 ) , FY ( 500 ) , FZ ( 500 ) , PAC , P FP TC , FM 
COMMON/C 0M7/ PPTC ,TPOT, PPE{ 1000) ,TLP E, ROEL , ROEL2, NEW 
C0MM0N/C0M8/R0EA , RO EB , ROEC , RO EC 2 , CP 0 ,CP 1 , C P2 ,C P3 , 

1 CFO , CF 1 ,CF2»CGD1 ,CGD2» CGB1, CGB2,CGF1,CGF2 
COMM CN/C 0M9/XN VAC, YNVAC, ZN VAC 
COMMON/ C OM 1 0/ I XN EW , IYNEW, IZNEW, II 

COMMON/C 0M1 1/RXNEW I (25 0 ) , RYN EW I { 2 50 ) , RZNEW I ( 250 ) , 
1KEEP (250) , NNUM( 250) 

COMMON/COMA/ A(4,5),MCR0 



LIST FORMAT STATEMENTS OF ALL READ COMMANDS. 



9010 FORMAT (2 0A4 ) 

9020 FORMAT ( 8A4, 3F8 . 5 » 2F 5. 2 ) 

9030 FORMAT ( 4 A4 , 3 F8 . 5 , 6 A4, F6 . 2 ) 

9040 FORMAT (F6.3*5X, I5,5I4,4X,3F5,3I2) 

9041 FORMAT (2 C 6.4,6F6.3) 

9042 FORMAT (F8.4,F8 .4 ) 

90 50 FORMAT ( 1 0A4 , A4 , 41 3 , F 8. 4 , 1 4 , F5 . 4 ) 
9052 FORMAT (6(F6.3)) 



LIST FORMAT STATEMENTS OF ALL WRITE COMMANDS. 



9610 FORMAT ( 1 HI 1 

9620 FORMAT ( 47X, ' SUMMARY OF A TOMS • / / , 3 5 X , 8A4 , • , NT ='I4 ,//, 
1 3 ( ' ATOM POSITION BIND ENERGY •),//) 

9630 FORMAT (3(1 5,3F6.2,F8.4,8X)) 

9640 FORMAT ( / 4X , F 10 • 3 , 2 5H EV, TOTAL KINETIC E N ERGY , , F 1 0 . 3 , 

1 27 H EV, TOTAL POTENTIAL ENERGY , F 10 . 3 , ' EV, REDUCTION',/ 
1/60X, 'RADIUS =' , F 5. 2 , i 

9650 FORMAT ( 105X , 4HP AGE , 13, /, 1H 1 ) 

9660 FORMAT ( / • ATOM DX DY DZ 

1VX VY . VZ KE PE TE',/) 

9670 FORMAT ( 118, 3F10.3,3F10. 1,3F10.4 ) 

9680 FORMAT ( ' SHARP DT D ECRE AS E • , 2 El 0 .3 ) 

9690 FORMATt 14, 3F5.2, 14) 

9691 FORMAT (9F8 .4 ) 

9692 FORMATt IX, 14,/) 

9693 FORMAT (4( I5,3X,F8.4,9X) ) 

9694 FORMAT ( 22X ,' SUMMARY OF POSITIVE POTENTIAL ENERGY CHAN 

1GES GREATER THAN',F7.4,2X,'NT=' ,12,//,' ATOM INITIA 
1L/FINAL X INITIAL/FINAL Y I N I T I AL /F I NA L Z PE 

1CHANGE ' /) 

9695 FORMAT ( 22X, ' SUMMARY OF NEGATIVE POTENTIAL ENERGY CHAN 
1 GES GREATER THAN • , F7 . 4 , 2X , • NT= ' , I 2 , // , • ATOM INITI 
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1AL/FINAL X INITIAL/FINAL Y I NI T I AL/ F I NAL Z PE 
1 CHANGE'/) 

9696 FORMAT ( 1 5 , 7X , 3 ( F6 . 2 , 1 3X ) , / , 12X , 2 ( F 6 . 2 , 13 X ) , F6. 2 , 
110X,F8.4) 

9697 FORMAT ( 10X , • I N I T I A L SUMMARY OF POTENTIAL ENERGY CHANG 

1ES ATOM PE CHANGE ')) 

9699 FORMAT (40X, 'FORCES AFTER EACH STEP',//,' STEP',10X, 
1'DT/STEP • , 11X , • FMAX* , 9X, • ATOM WITH MAX FORCE', 5X, 

1 ' FORCE (1)'/) 

9700 FORMAT ( 2X, I 3 , 6X , E 1 2 . 4 ,6X , E 12 . 4 , 1 3X , I 3 , 12 X , E 1 2 . 4 ) 

9704 FORMAT (20X,'DT CHECK FOR EACH STEP ',//,' STEP ', 1 IX , 

1 ' STEP ' ,1IX,'DTE1' *11X, 'DTF1' ,11X,'DTE',11X,'DTF',//J 
9704 FORMAT ( 20X , ' DT CHECK FOR EACH STE P ',//,' ST E P ', 1 1 X , 

1 ' DTE1' ,11X, 'DTF1 • , 11X, 'DTE • , 11X, 'DTF* ,// ) 

jJ: ){c % j(r * £ £ A £ * :$( # £ 'Jr- ^ ^ ^ ♦ 3(c s}: s{i & Jfc # # # 

C 'FORCE CONFIGURATION' 

J$C # if. if. # ^ # jfc * # && £ ;}>: 3^ ^ # ;£ # rfs £ £ * # # # £ 



CHANGE STATEMENT NUMBER 9704 TO READ 



9704 FORMAT ( 20X, ' DT CHECK FOR EACH STEP ',//• STEP ', 1 IX , 
1 'DTF ' , 12X, 'DTI A' ,11X, 'DTE' »11X, ' DTIV' ,/ / ) 

9705 F0RMAT(I5,4(3X,EI2.4) ,/) 



INITIALIZE APPROPRIATE VARIABLES AND VARIABLE ARRAYS, 



ST ART = 0 . 0 1 * I T I M E ( XX ) 

DO 2 1=1,1000 

RXK ( I ) = 0 . 0 

RYK ( I ) = 0 .0 

RZK( 1 ) =0.0 

VX( I ) = 0. 0 

VY ( I ) = 0 .0 

VZ( I ) =0. 0 

PK E ( I ) = 0 .0 

ppe ( n=o.o 

PTE ( I ) =0, 0 
2 RZ I ( I ) = 0 . 0 
I SHUT = 1 
NRUN= 0 



READ INPUT DATA COMMON TO ALL DESIRED CALCULATIONS. 



READ ( 5,9010) 

READ ( 5,9020) 
READ ( 5,9030) 

READ ( 5,9030) 
READ ( 5,9050) 

READ (5,9042) P 



I H 1 

I H2 , DCON, ALPHA, RE,ROEC, ROEL 
BULLET, 8 MAS,PEXA,PEXB,IHB, THE RM 
TARGET ,TMAS,EXA,EXB, IHT, TEMP 
IHS, PLANE, LS, IX ,IY, IZ,CVR,MCRO , 
P E PC K , P P ENC K 



DTI 



INITIALIZE CONSTANTS AND DEFINE SCALING FACTORS. 
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DT I S = DT I 
R0E2=4. 0 
R0E=S0RT(R0E2) 

ROEM = ROE-DTI 

ROEL2 = ROEL*ROEL 

CVE= 1. 60E-19 

CVM=1.672E— 27 

VF AC =0. 5 

FM= 1 . 0E- 10 

FM2 = FM^FM 

C VO =CVR * 1 . 3E-1 0 

CVED=CVE/CVD 

PTMAS=TMAS*CVM 

PBMAS=BMAS*CVM 

HTMAS=0.5^PTUAS/CVE 

HBMAS=0 .5*PBMAS/CVE 

TSAVE=BMAS/( BMAS+TMAS) 

BSAV E=TMAS/ ( BMAS + TMAS) 



DEFINE REPULSIVE POTENTIAL PARAMETERS. 



FXA = AL0G t-EXB« CVED) +EXA 
PFXA=ALQG(-PEXB*CVED)+PEXA 
PPTC=EXP (PEXA+PEXB*ROE ) 
PAC=ALOG(CVED) +PEXA 
PFPTC = EXPt PAC + PEXB*ROE ) 



DEFINE ATTRACTIVE POTENTIAL PARAMETERS. 



CGD1=AL0G( DCOM ) +2.0* ALPHA*RE 

CG02=AL0G(2.D*DCCN)+ALPHA*RE 

CGBl=-2. 0*ALPHA«CVR 

CGB2=— ALPHA^CV R 

CGF 1 =ALOG (-CGBl^CVED) +CGD1 

CGF2 = AL0G( -CGB 2*CVE D) +CGD2 



DEFINE REGIONS OF APPLICABILITY OF POTENTIAL FUNCTIONS 



ROE A=1 . 50/CVR 
R0EB = 2. 0/C VR 
R0EC2=R0EC*R0EC 



DEFINE PARAMETERS USED TO DETERMINE THE BEST CUBIC 
FIT BETWEEN THE MAXIMUM DISTANCE OF APPLICABILITY OF 
OF THE REPULSIVE POTENTIAL AND THE MINIMUM DISTANCE 
OF APPLICABILITY OF THE ATTRACTIVE POTENTIAL. 
SUBROUTINE CROSYM PERFORMS THE NECESSARY CALCULATIONS. 



A( 1 1 1 ) = 1 .0 

A { 1 1 2 ) =ROE A 

At 1 ,3) =ROEA*ROEA 

At 1,4) = RCEA**3 

All ,5) =EXP( EX A+ EX B*ROEA) 

At 2. 1 ) = 1 . 0 
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A( 2 * 2 ) = ROE B 
A ( 2 » 3 ) = ROE B A ROEB 
A( 2,4) =R0EB**3 

A (2 , 5)=EXP (CGD1+CGB1*R0EB)-EXP( CGD 2+CGB 2*R0E B ) 

A(3,l) =0.0 

A( 3, 2) =- 1.0 

A(3, 3)=-2 .0*R0EA 

A( 3,4) =— 3 . 0* ROE A* ROE A 

A( 3, 5)=EXP(FXA+EX6*R0EA) /CVED 

A ( 4 , 1 ) =0 .0 

A ( 4 , 2 ) =-1.0 

A (4, 3) =-2. O'^ROEB 

A ( 4 , 4 ) =— 3 .0*RGE3*R0EB 

A( 4,5) =( EXP (CGF1+CGB1* ROEB) -EXP (CGF2+CGB2* ROEB) )/CVED 
CALL CROSYM 
C P0 = A ( 1 ,5) 

CP 1= A ( 2 , 5) 

CP 2= A ( 3 , 5) 

CP3=A (4,5) 

CF 0=-CP 1*CVED 
CF1=-2.0^CP2*CVED 
CF2=-3. 0*CP3*C VED 



READ INPUT DATA FOR EACH SITE TO BE INVESTIGATED. 
MULTIPLE INVESTIGATIONS ARE POSSIBLE BY SIMPLY READING 
IN DATA FOR MORE THAN ONE SITE. 



5 READ ( 5,9040) EVR , NTT , NS , ND , I P , I DE EP , I T Y P E , D1X, 
1D1Y* D1Z, I VAC X , I VAC Y , I VACZ 
RE AD (5, 9041 ,END=99 99) DDT I I , DDT I F , XNVAC , YN VAC , ZN VAC 

5?C^C 34c si:## * 3$C;$C £ If * * ^3>C * # jfif A * 5*t # # ** ## # £ # 3$C # * 5& # 2$C # # £ 3$C jfc 

•FORCE CONFIGURATION' 



IN THE FORCE CONFIGURATION ONLY, READ IN DTI VALUES 
TO BE ASSIGNED AFTER FORCE COMPARISON, AND 
DEFINE VARIABLE VFAC2. 



READ (5,9052) 0TIA1,DTIA2,DTIA3,DTIV1,DTIV2,DTIV3 
VF AC2 = VFAC*V FAC 

# 3}C 3*C5£ # # £ # * # jfc £ # # # * * & # 5$C # # £ # £ >jC £ £ * £ £ £ £ £ 



IF(NTT.EO.O) GO TO 9999 
IQ= ITYPE— 1 
EV=EVR*1.0E+3 
DT I = DT I S 
TPKE =EV 



sjcsjt 3c £ £ sfcAafcsScc'c ££ jjcijcjjcA^csjcsic 3)e 3: cfccfcijcijc jje#* 3(t*4c^£*^:^:j'{3)e3jc^c^^:^:^c^c*tX c:, !'^ : * :s P* 



sjc^c^eAcfcioJofc^ujccfccfcjjc^tc^^ejkstcjJc^c^c-jc-bjjtsjcciesfe^si-rfcs^jSc^c^cije 

•POTENTIAL WELL CONFIGURATION’ 

Jjojojc^cstofc^cje^cilcjfc^ccijcjjc* jjc^csjejjc^r :fcj{c:i(;:>5c5io)e3c*3«*'fc 



CONSTRUCT A 'RELAXED’ CRYSTAL IN THE COMPUTER BY 
READING IN ’RELAXED* CRYSTAL PARAMETERS AND POSITIONS 
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VICE USING SUBROUTINE B100 TO CONSTRUCT THE CRYSTAL. 



READ (5,9690) LL , D1X , D 1Y , D 1Z , NV AC 
DO 15 I =1 , L L , 3 
K= I + 1 
J=I+2 

READ (5,9691) RX ( I ) , RY ( I ) , RZ ( I ) , RX ( K ) , RY ( K ) , RZ ( K ) , 
1 RX ( J ) , RY ( J ) , RZ ( J ) 

15 CONTINUE 



•OTHER CONFIGURATIONS' 



SELECT THE DESIRED CRYSTAL STRUCTURE AND ORIENTATION. 
SUBROUTINE B100 CONSTRUCTS THE (100) PLANES OF A BODY- 
CENTERED CUBIC CRYSTAL IN THE COMPUTER. 



14 CALL B100 

** * *** * *************************** *********************** 



30 ILAY=IDEEP 

IF(IDEEP) 35,35,40 
35 LD=LL 
I L AY = I Y 
40 RLL=1. O/LL 
TPOT L= 1 .0 
DO 45 I =1 , LL 
RXK ( I ) = RX( I ) 

RYK ( I ) = RY ( I ) 

RZK ( I ) =R Z ( I ) 

RXI ( I ) =R X{ I ) 

RYK I ) = RY ( I ) 

45 RZI ( I ) =RZ ( I J 

IF(NRUN.EO.O) GO TO 60 
DO 55 I =1 , LL 
LCUT ( I )=0 
RX ( I ) = PX Id) 

RY ( I ) =RYI { I ) 

RZ ( I ) = R Z I ( I ) 

RXK ( I ) = RX I ( I ) 

R YK{ I ) =R YI ( I ) 

55 RZK ( I ) = RZI ( I ) 

60 NRUN=1 



THIS SECTION CALCULATES THE ENERGIES OF ALL ATOMS IN 
THEIR INITIAL POSITIONS IN THE PERFECT LATTICE (THAT 
IS, WITH NO INTERSTITIAL IMPLANTED). INITIAL POSITIONS 
AND ENERGIES OF ALL ATOMS ARE PRINTED TO PROVIDE A 
COMPARISON WITH CHANGES IN CRYSTAL ATOM POSITIONS 
AND ENERGIES CAUSED BY IMPLANTATION OF THE INTERSTI- 
TIAL. 



************** ********************************************** 



C 



•POTENTIAL WELL CONFIGURATION' 
************ ************ *********** 



c 

C THIS CALCULATION CAN NOT BE MADE IN THE POTENTIAL 

C WELL CONFIGURATION WITHOUT DESTROYING THE INPUT 

C DATA OF THE 'RELAXED' CRYSTAL. CONSEQUENTLY, 

C 

C R X( 1 ) = 25 . 0 

C 

C SHOULD BE DELETED IN THE POTENTIAL WELL CONFIGURATION. 

C 

;3c # # # rfr # % # # £ $ £ # # & £ # £ # # £ £ # # # £ Of. £ # if. D$C 3$G # # 5$C # £ # 3$C # # # afc # # # # # 3jc 3$C 35c # 5ic 



RX ( 1 ) =2 5 . 0 
DO 63 1=1, LL 
VX ( I )=0.0 
VY( I ) =0. 0 
VZ( I )=0.0 
PPEl I )=0 .0 
PKE ( I ) =0.0 
63 PT E ( I ) = 0 .0 
TPOT =0.0 
N E W = 0 

CALL ENERGY 
NPAGE=1 
NT = 0 

WRITE ( 6,9610) 

WRITE ( 6,9620) IH2,NT 
DO 61 1=1, LL, 3 
K= I +1 
J = I + 2 

61 WRITE ( 6,9630) I , RX( I ) , RY ( I ) , RZ (I) , PPE (I) , K ,RX ( K) , 
1RY(K),RZ(K),PPE(K),J,RX(J),RY(J),RZ(J),PPE(J) 

WRITE ( 6,9650) NPAGE 
DO 62 1=1, LL, 3 
K= I +1 
J = I +2 

PPFS AV ( I ) = PPF( I ) 

PPES AV { K) =PPE( K) 

PPE S AV ( J ) =P PE ( J) 

62 CONTINUE 



3$C5$r # ## # # £ * sJl # 3*C 5$C 

C 'POTENTIAL WELL CONFIGURATION' 

58c 5§C 5?C 5}C 5^ 5}C 5^ 



c 

C THE VARIABLES BELOW ARE USED TO CREATE OFFSETS FROM 

C EQUILIBRIUM POSITIONS IN THE 'RELAXED' CRYSTAL. IF 

C NO OFFSET IS DESIRED, THESE VARIABLES SHOULD BE 

C INCLUDED BUT SET EQUAL TO ZERO. THE OFFSET IS 

C INSTITUTED BY SUBROUTINE PLACE. THESE VARIABLES 

C SHOULD NOT BE INCLUDED IN GTHER CONFIGURATIONS. 

C 



D1X = 0 .0 
D1Y=D. 0 
D1Z= 0 . 0 



c 

C SUBROUTINE PLACE CREATES THE DESIRED VACANCY, INTER- 

C STITIAL, OR SELF INTERSTITIAL IN THE CRYSTAL. 

C IN THE POTENTIAL WELL CONFIGURATION, THE INTERSTITIAL 
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IS OFFSET FROM ITS EQUILIBRIUM POSITION. 



CALL PLACE 
DO 6 5 I = 1 , L L 
VX ( I ) =0 .0 
VY ( I ) =0. 0 
VZ ( I ) =0 .0 
PPE ( I ) =0 .0 
PKE ( I ) =0. 0 
65 PT E ( I ) = 0 . 0 
TP0T=0 . 0 
NE W = 0 



THE ENERGY SUBROUTINE NOW CALCULATES ENERGIES OF ALL 
ATOMS OF THE CRYSTAL AFTER IMPLANTATION. THESE 
ENERGIES AND THE INITIAL POSITIONS OF ALL ATOMS ARE 
PRINTED FOR TIME ZERO. CHANGES IN POTENTIAL ENERGY OF 
ALL ATOMS AS A RESULT OF IMPLANTATION ARE ALSO CALCU- 
LATED AND PRINTED. 



CALL ENERGY 
B I ND=— T POT 
TE=TPOT+BI ND 

T I ME=0 . 0 
NT =0 

WRITE ( 6.9620) IH2.NT 
DO 70 1=1, LL, 3 
K = I + 1 
J = I + 2 

70 WRITE ( 6,9630) I,RX(I),RY( I ) , R Z ( I) ,PP E ( I ) ,K,RX(K) , 
1RY(K),RZ(K),PPE(K),J,RX(J),RY(J),RZ(J),PPE(J) 

WRITE ( 6,9640) TPKE , TPOT , T E , ROE L 
NPAGE= NP AGE + 1 
WRITE ( 6,9650) NPAGE 
WRITE (6,9697) 

DO 80 1=1, LL, 4 
K = 1 + 1 
J= I + 2 
L= 1+3 

PPE INT ( I ) =PPE( I)-PPESAV(I) 

PPEINTt K ) = PPE( K )-PPESAV( K) 

PPE I NT ( J ) = PPE( J )-PPESAV( J ) 

PPEINT( L)=PPE( L)-PPESAV(L) 

80 WRITE (6,9693) I , PP E I NT ( I ) , K , P PE I NT ( K ) , J , PP E I NT ( J ) , L , 
1 PPE I NT ( L ) 

NPAGE =NPAGE+ 1 
WRITE (6,9650) NPAGE 
DT =1 . 0 E-l 5 



•FORCE CONFIGURATION* 



SINCE TIMESTEPS IN THIS CONFIGURATION ARE DETERMINED 
DIRECTLY BY A FORCE COMPARISON, NO TIMESTEP 
DECREMENTING PROCEDURE IS REQUIRED. THE ONLY CARD 
NEEDED IS 



100 CONTINUE 
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c 



sic###*## ^c^c5k>jc^r:>}c*3k3jc*3^5jc5t:5?c^c^c5?:^^^c5!c^:^c^:3*c>jc^c5!f 

•OTHER CONFIGURATIONS' 



THE FOLLOWING STATEMENTS DEFINE THE DTI 
DECREMENTING PROCEDURE. 



DDTI =DDT I I 
NDEC =0 

100 DTI = DTI-DDTI 

IFIDTI. LT.DDTIF) DT 1= DDT I F 
105 NDEC=NDEC+ 1 



THE MAIN BODY OF THE PROGRAM NOW SOLVES THE EQUATIONS 
OF MOTION BY THE AVERAGE FORCE METHOD AND DETERMINES 
POSITIONS OF ALL ATOMS AT THE END OF THIS TIMESTEP. 
SUBROUTINE STEP PERFORMS ALL FORCE CALCULATIONS. 



DTOD=DT /CVD 

TFAC=2.0*PTMAS*DTI *CVD 
TFAC8=2. 0*PBMAS*DT I*CVD 
TEFAC=DTI*CVD 
HDT0D=0 ,5*DT0D 
DTOM = DT / PTMA S 
HDTOM= 0 . 5*D TOM 
DTOMB = DT/P BMAS 
HD TO MB = 0. 5* D TOMB 
200 CALL STEP 

IF (LCUT(l) .GT.O) GO TO 240 
1 = 1 

RXK ( I ) = RX ( I ) 

RYK ( I ) =RY{ I ) 

RZK ( I ) = RZ ( I ) 

RX ( I ) = RX ( I ) +DT OD* (HDTOMB*FX( I ) + VX ( I ) ) 

RY( I ) =RY( I )+DTCP*(HDTOMB*FY ( I )+VY ( I ) ) 

RZ( I ) =RZ ( I )+DTOD*(H!)TOMB*FZ< I ) + VZ ( I ) ) 

240 DO 245 I =2 , L D 

IF( LCUT( I ) .GT.OJGO TO 245 
RXK ( I ) = R X ( I ) 

RYK ( I ) =RY ( I ) 

RZK( I ) =RZ ( I ) 

RXC I )=RX( I )+DTOD*(HDTOM*FX( I )+VX( I ) ) 

RYU )=RY(I ) +DTOD* ( HDTOM*FY ( I )+VY( I ) ) 

RZ ( I ) =P Z ( I I +DT OD* ( HDTOM*FZ( 1 J+VZ ( I ) ) 

245 CONTINUE 
CALL STEP 

FSTACC(NT)=SQRT(FX( 1)**2+FY( 1 ) ** 2+F Z ( 1 ) ** 2 ) / BMAS 
EMAX = 0 .0 
FMA X = 0. 0 
T IME = T I M E+DT 
NT =NT + 1 

IF<LCUT( 1) .GT.O) GO TO 265 
1=1 

VSS=VX ( I ) 

VX( I ) = VSS+HDTOMB*FX (I ) 

RX( I )=RXK ( I )+( VX( I ) +V S S ) *HDTOD 
VSS = VY ( I ) 

VY ( I ) =VSS+HDTOMB*FY (I ) 

RY( I ) = PYK ( I ) + { VY ( I )+VSS)*HDTOD 
VSS=VZ ( I ) 



73 



oonoo 



VZ ( I ) =VSS+HDTOMB*FZ(I ) 

RZ( I ) = P. Z K ( I ) + { V Z ( I )+VSS)*HDTOD 

PKE( I )=VX( I )*VX( I ) +VY ( I ) * VY ( I ) + VZ ( I ) *VZ ( I ) 

EMA X 1 =P KE { I ) 

FMAX 1 2= FX ( I )*FX{ I)+FY( I )*FY( I ) + FZ(I )*FZ( I ) 
FMAX1 = SQRT ( FMAX12 ) 

AMA Xl= FMAX 1 /BM A S 
FSSACC ( NT ) = AMAX1 
260 F X ( I ) =0.0 
FY ( I 1 = 0.0 
FZ ( I ) = 0 .0 
FMAX=0. 0 
EMAX= 0.0 
F2 M=0 . 0 

265 DP 280 I =2, LD 

IF(LCUT(I).GT.O)GO TO 280 
VSS=VX( I ) 

VX{ I ) =VSS+HDTOM*FX(I } 

RX ( I ) = RXK ( I ) + ( VX ( I ) +VSS J^HDTOD 
VSS=VY ( I ) 

VY ( I ) = V S S+HUTO M*F Y ( I ) 

RY ( I ) = RYK { I ) + ( VY ( I J+VSS )*HDTOD 
VSS=VZ( I ) 

VZ( I ) = VS S + HDTD M^-'F Z ( I ) 

RZ ( I) = RZK ( I ) + (VZ ( I ) +VSS )*HDTOD 
PKF ( I)=VX(I)*VX(I)+VY(I)*VY(I)+VZ(I)*VZ(I) 
275 F2 = FX( I )*FX< I)+FY( I )*FY< I)+FZ( I )*FZ( I ) 

FX ( I ) = 0. 0 
FY ( 1 ) = 0 . 0 
FZ ( I )=0 .0 

I F ( F2. LE.F2M) GO TO 278 
F2M=F2 

I FMX AM ( NT ) = I 

278 IF ( PKE ( I } . GT. EMAX) EMAX=PKE(I) 

280 CONTINUE 

FMAX=SQRT(F2M) 

AMAXL=FMAX/TMA S 
FORMAX { NT ) = AMAXL 



C 'FORCE CONFIGURATION' 



TIMESTEP DETERMINATION IN THE FORCE CONFIGURATION 
IS PERFORMED BY COMPARING THE MAXIMUM FORCE IN THE 
CRYSTAL WITH APPROPRIATELY CHOSEN TEST VALUES. 



AMAX= AMAX1 

IF(AMAXL.GT.AMAX) GO TO 282 
I FMXAM ( NT ) = 1 
GO TO 286 
282 AMA X=A MA XL 
2 86 DTL = DT 

EMAXL = E MAX^V FAC2 

IF( AMAX. LE. l.OE-8) DTI A = DTI A1 

IF(AMAX.LE .1.0E-9) DTI A=DT IA2 

IF (AMAX.LE.1.0E-10) DT I A= DT I A3 

I F ( AMAX.LE. l.OE-ll) ISHUT=-1 

CT IME=0.01*ITIME<XX)-START 

EMAX=EMAX1 

I F ( EMAXL .GT.EMAX) EMAX = EMAXL 
I F( EMAX .LE.l .0E+6) DTIV=DTIV1 
I F { EMAX. LE. l.OE+6) DTIV=DTIV2 
IFtEMAX.LE. l.OE+2) DTIV=DTIV3 
DTFCK= ( DTI A*CV D*2 .0 )/ ( AMAX/CVM ) 
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IF(DTFCK.LT.C. 0) DTF CK =-DTFCK 
DT F=S0RT (DTFCK ) 

DTF S P ( NT ) = DT F 

DTI ASP( NT) =DTI A 

T F { EMAX .GT .0.0 ) GO TO 290 

DTE=1 . OE-5 

GO TO 295 

290 DTE=(DTIV*CVD)/SQRT(EMAX) 

295 DTE SP( NT ) =DTE 

DTIVSP(NT)=DTIV 

DT=DTF 

I F ( DTE . LE.DT)DT=DTE 
DTSTEP(NT)=DT 



C ’OTHER CONFIGURATIONS' 

£s!csjc£:$-fc £ A s£ jfr. aJcjjcjfc Jf- :}= A 



FOUR NEW TIMESTEPS ARE CALCULATED BASED ON MAXIMUM 
FORCES AND ENERGIES OF INITIAL AND FINAL POSITIONS. 
THE SMALLEST IS CHOSEN AS THE NEXT TIMESTEP INTERVAL. 



DTL=DT 

CTIME=0. 01*1 TIME ( XX) -START 
DTE1 =TE F AC*SQRT ( 1 .O/EMAXI) 
DTF1=SCRT( TFACB/FMAX1 ) 
DTE=TEFAC*SQRT ( I .0/ EMAX) 

DTF = SQRT (TFAC/ FM AX ) 

DTI ( NT ) =DTE 1 
DT2 ( NT ) = DT F 1 
DT3 { NT ) =DT E 
DT4( NT) =DTF 

I FI EMAX1 .GT.EMAX) EMAX=EMAX1 
DT =DTE I 

IF(DT.GT.DTFl) DT=DTF 1 
IF(DT .GT .DTE ) DT = DT E 
IF(DT.GT.DTF) DT=DTF 
DT STEP { NT ) =DT 



DAMPING IS INTRODUCED IN THE FORM OF A DAMPING FACTOR 
WHICH DECREMENTS ALL VELOCITY COMPONENTS. 



300 I F { I S HUT . E 0 . — 1 ) GO TO 400 
310 IF(NS-NT) 400,400,320 
320 DO 325 1=1, LL 

VX( I )=VFAC*VX< I ) 

VY ( I ) = VF AC* VY ( I ) 

325 VZ ( I ) =VFAC*VZ( I ) 



% nje :$c sic 3p £ sicsic # s& sic afc $ * sic * sje# >^5^^c3jc){c^'>!;^3^3t^C3}c^;3{r^:3fccxt3}c^^ac3!crfc5^3}:^3jt3t:^:^5^A^ 

’FORCE CONFIGURATION’ 



SINCE THE TIMESTEP HAS ALREADY BEEN DETERMINED IN THE 
FORCE METHOD, THE COMPUTATION CAM NOW SHIFT TO 
STATEMENT NUMBER 100 AND BEGIN THE DYNAMICS FOR THE 
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NEXT TIMESTEP. 



GO TO 100 



Jfc # # # # # a?c # 3$C * * 3?C 3*C 5?C 2{C 5»c & A 3*05* 3$C 3$C # * 5*C # * ^ 2C # # 5$C 

'OTHER CONFIGURATIONS' 

jfr 3 {c # * jjc Jff A sjs J)c # 5$c ^ ^ + ♦ * >)c s!t A ^ ❖ * # £ * £ % 



SHIFT TO STATEMENT NUMBER 800 FOR CALCULATION OF THE 
DTI DECREMENT IF REQUIRED FOR THIS TIMESTEP. 

GO TO 800 



!$C 5k if. J$C 5}c # 3*C if. 5fc * 3* # 3$< # 3$C £ 3{c 3fc # ^ 



:fe#5£3}c 5}c* jJcsJc# j5c:fc 5fe^c^c^c5}t3jc5}c^c>{c^c^:3{c5?c^c^C3jc^ok^: 



SUBROUTINE PRINT PRINTS PERTINENT GENERAL DATA FOR 
EACH TIMESTEP PRINTOUT. 



400 CALL PRINT 



RELATIVE MOTION, VELOCITY, AMD ENERGY OF EACH ATOM 
ARE PERIODICALLY PRINTED. 



410 TPOT = Q. 0 

DO 450 1=1, LL 
PP E ( I )=0 .0 
450 PTE( I ) =0. 0 
CALL ENERGY 
PKE ( 1 ) = HBMAS*PKE( 1 ) 

TPKE =PKE ( 1 ) 

PT E ( 1 ) = PKE ( 1 ) +PP E ( 1) 

DO 620 I =2 , L L 
P KE ( I ) =HTMAS*PKE ( I ) 

TPKE=T PKE+PKE ( I ) 

620 PTE ( I ) = PKE ( I ) +PPE( I ) 

TE=TPOT+BI NO 
WRITE ( 6,9660) 

DTEST= 0. 1=M RX ( 1 ) -RX I C 1 ) )**2 
IF (DTEST.GT. 0.01) DTE ST= 0.01 
IF(TPOT.LE.TPOTL) GO TO 700 
ER AT =TP KE*RLL 
7 00 DO 7 50 I = 1 , LD 

DX ( I ) = RX ( I ) -RX 1(1) 

DY ( I ) =R Y ( I ) — RY I ( I ) 

DZ( I ) = RZ ( I )-RZ 1(1) 

IF ( DX ( I )*#2.GE.DTEST ) GO TO 720 
IF ( DY( I ) * ~ 2 . G E .DTEST) GO TO 720 
IF (DZ( I )**2.GE. DTEST ) GO TO 72C 
GO TO 750 

720 WRITE ( 6,9670) I , DX ( I ) , DY ( I ) , DZ ( I ) » VX ( I ) , VY( I ) , 
1VZ( I ) , P KE ( I ) ,PPE ( I ) , PTE (I ) 

750 CONTINUE 

WRITE ( 6,9640) T PKE , T POT , T E, RO EL 
NP AGE=NP AGE+ 1 
WRITE ( 6,9650) NP AGE 
TPOT L=T POT 
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IF(NT-NTT) 760,950,950 
760 DO 780 1=1, LL 

VX ( I ) =V F AC*V X ( I ) 

VY ( I ) =VFAC~VY( I ) 

VZ( I ) = V P AC*V Z ( I ) 

780 CONTINUE 

IF( I SHUT.E0.-1) GO TO 950 
790 NS=NS+ND 



ijc skjfe # 5jc # # jj: ^ 5jc if. # * # * # + s)c # # # + if if ♦ * ^ ;jt A :£ # jjc # £ j{; it ^ # # * # # ^ jje 

C . 'FORCE CONFIGURATION' 

3^ # if. # ajs jfr: # # # * 5*; •* :£ # # # 5*; # * 3j< * * *. # £ # # * 5}c 



NO DECREMENTING PROCEDURE FOR DTI IS NEEDED, SC THE 
FORCE CONFIGURATION SHIFTS TO STATEMENT NUMBER 100 
AND BEGINS COMPUTATIONS FOR THE NEXT TIMESTEP. 



830 GO TO 100 



^ * & 2$c :& # 

•OTHER CONFIGURATIONS' 



THE DTI ALTERATION PROCESS IS BEGUN FOR THE NEXT 
TIMESTEP 



800 IF (NDEC. EO.l 0) GO TO 810 
GO TO 820 
810 DDT 1 = 0.1 *DDT I 
DT I =DT I + DDT I 
N 0 E C = 0 

820 GO TO 100 



950 CONTINUE 



FINAL POSITIONS (IN RECTANGULAR COORDINATES) AND 
BINDING ENERGIES OF ALL ATOMS ARE PRINTED AFTER THE 
LAST TIMESTEP. WRITE (7,XXXX) STATEMENTS ARE INCLUDED 
WHEN DATA DEC*S CONTAINING COMPLETE INFORMATION FOR 
THE ENTIRE CRYSTAL ARE DESIRED. ADDITIONALLY, POSITIVE 
AND NEGATIVE POTENTIAL ENERGY CHANGES GREATER THAN 
A PRESET VALUE ARE DETERMINED AND PRINTED. 



955 WRITE ( 6,9620) IH2,NT 

WRITE (7,9690) L L , D1X , D1Y , D1 Z , NV AC 
DO 965 1=1, LL, 3 
K= I +1 
J = I +2 

WRITE ( 7,9691) RX( I ) , RY( I ) , RZ ( I ) , RX ( K) , RY ( K) , RZ ( K ) , 
1RX(J) ,RY( J),RZ( J) 

965 WRITE ( 6,9630) I , RX ( I ) , RY ( I ) , RZ ( I ) , PP E { I ) , K , RX ( K ) , 
1RY(K) , R Z ( K ) ,PPE(K) ,J,RX(J) ,RY(J) ,RZ( J) ,PPE(J) 

WRITE ( 6,9640) T PKE , TPOT , T E , RO EL 
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NPAGE =NP AG E+ 1 
WRITE ( 6,9650) NPAGE 
DO 970 I =1 , LL 
PPEPOS ( I ) = 0. 0 
PPENEGC I ) = 0.0 

PPKEEPII )=PPE( I)— PPESAV(I) 

IF(PPKEEP( I ) .ST.PPEPCK) GO TO 968 
IF(PPKEEP( I ).3T.PPENCK) GO TO 970 
PPENEG ( I ) = P PKE E P ( I ) 

GO TO 970 

968 PPEPOS ( I ) = P PK E EP (I ) 

970 CONTINUE 

WRITE (6,9694) PPEPCK.NT 
DO 980 1=1 ,LL 

IF(PPEPOS( I ) .LT. PPEPCK) GO TO 980 

WRITE ( 6,9696) I , R X I ( I ) , RYI ( I ) , RZ I ( I ) ,RX ( I ) , RY( I ) , 
1 RZ( I ) , PPEPOS ( I ) 

980 CONTINUE 

NPAGE=NP AGE+1 
WRITE (6,9650) NPAGE 
WRITE(6,9695)PPENCK,NT 
DO 990 1=1, LL 

IF ( PPENEGl I ) .GT. PPENCK ) GO TO 990 
WRITE ( 6,9696) I , R XI ( I ) , RYI ( I ) , RZI ( I ) , RX ( I ) , RY ( I ) 
1RZ( I ) ,PPENEG( I ) 

990 CONTINUE 

NPAGE=NPAGE+1 
WRITE (6,9650) NPAGE 
WRITE (6,9699) 

DO 995 1=1. NT 

995 WRITE (6,9700) I , DTST EP ( I ) , FORM AX ( I ) , I FM XAM ( I ) , 

1 FS SACC ( I ) 

NPAGE=NP AGE+1 
WRITE (6,9650) NPAGE 
WRITE (6,9704) 

DO 999 1=1, NT 






C 'FORCE CONFIGURATION' 



DATA PERTINENT TO TIMESTEP DETERMINATION BY THE 
FORCE METHOD ARE PRINTED. 



999 WRITE ( 6,9705) I ,DTFSP ( I ) , DTI A SP ( I ) ,DTESP ( I ) , DT I VSP ( I ) 



ak :fc afc * afc afc # ak # 

•OTHER CONFIGURATIONS' 

afc###a^£>{ca!c3!saic2$ca$e afraic# a^akajc# £ a?:## ajcafcafc a(c # a$c £ * aAc 



DATA PERTINENT TO "TIMESTEP DETERMINATION BY THE 
ENERGY METHOD ARE PRINTED. 



999 WRITE (6,9705) I , DT 1 { I ) , DT 2 ( I ) , DT 3( I ) , DT4 ( I ) 
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c 



• PLUCK CONFIGURATION' 

# 3jc^C^5}Cil:^C 3^# *^>{C^^^lC5JC3k3}S>{C^:>^5k>^>JC3*C5jc3?S3i:3^^:>{c3fe3^5^^ 



SUBROUTINE PLUCK CAN BE USED TO PUNCH DATA CARDS FOR 
A SMALLER CRYSTAL CENTERED ON THE INTERSTITIAL FOR 
USE IN THE DYNAMIC PROGRAM. DELETE THESE CARDS IF 
DATA FOR THE SMALLER CRYSTAL IS NOT DESIRED 



CALL PLUCK 

WRITE (7,9690) L L , D 1 X , D1 Y , D1 Z , NV AC , I XNE W , I YNE W , I ZNEW 
LL= I I 

DO 1100 1=1, LL 

1100 WRITE (7,9691) I , R XNE WI ( I ) , RYNE WI ( I ) , RZNE W I ( I ) , KE EP ( I ) 

jjcjfcifc jjcafc## sjc^ajc 



1000 I F ( ISHUT ) 9999,5,5 
9999 STOP 
END 



SUBROUTINE CROYSM SOLVES M SIMULTANEOUS EQUATIONS BY 
THE METHOD OF CROUT IN ORDER TO c I T THE BEST CUBIC 
EQUATION BETWEEN THE REPULSIVE AND ATTRACTIVE PARTS 
OF THE POTENTIAL. 



SUBROUTINE CROSYM 



COMMON/COMA/ A(4,5),MCR0 
M=MCRO 
N=M + 1 
11=1 

100 13=11 

SUM= ABS ( A ( 11,11)) 

DO 120 I =1 1 , M 

IF(SUM-ABS(A(I ,11) ) ) 110,120,120 

110 13=1 

SUM= ABS ( A ( I , II ) ) 

120 CONTINUE 

I F ( 13— II) 130,150,130 
130 DO 140 J=1 ,N 
SUM=-A ( I 1, J ) 

A ( 1 1 , J ) = A ( 13, J ) 

140 A(I3,J)=SUM 
150 13=11+1 

DO 160 1 = 13, M 

160 AC I ,11 ) =A( I ,11 )/A( I 1,1 1) 

170 J 2= 1 1- 1 
13 = 11+1 

I F ( J 2 ) 180,200, 180 

180 DO 190 J = I 3 , N 
DO 190 1=1, J2 

190 A( I 1 , J) =A( I 1 , J ) —A ( I 1 , I )*A( I , J) 
IF(Il-M) 200,220,200 
200 J2 =1 1 

I 1 = 1 1+1 
DO 210 I = I 1 , M 
DO 210 J=1 , J2 

210 A(I,I1)=A(I,I1)-A(I,J)*A(J,I1) 
IFII1-M) 100,170,100 
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. 220 DO 240 I =1 , M 
J2=M- I 
13 = J 2 + 1 

A(T3,N)=A( I 3 * N ) / A ( 13,13) 

I F ( J 2 ) 230,250,230 

230 DO 240 J = 1 , J2 

240 A( J,N) =A( J,N)-A( I3,N)*A(J,I3) 
250 RETURN 
END 



SUBROUTINE B1DD GENERATES A BODY-CENTERED CUBIC 
LATTICE IN THE (100) DIRECTION. 



C 



SUBROUTINE B100 



C0MM0N/CCM1/RX ( 500) ,RY ( 500 ) , RZ ( 500) , LCUT ( 500) , 
ILL ,LD, I TYPE ,NVAC 

COMMON /COM 4/ IX, I Y, I Z, SCX, SC Y , SC Z , I DE EP , D 1 X , D1 Y, 
1D1Z, I V ACX , I V ACY , IVACZ 
DIMENSION YLAX (20) 

DATA YLAX/ 20*0 . 0/ 

YLAX(l) =-0.20 
YLAX( 2) =-0. 03 
SCX = 1 .0 
SCY=1. 0 
SCZ=1.0 
M = 2 
JT =0 
Y=-SCY 

DO 60 J = l, IY 
Y = Y+ SC Y 
KT = 0 
Z=-SCZ 

DO 59 K=1 , I Z 
Z=Z+SCZ 
IT = 0 
X=-SCX 

DO 58 1 = 1, I X 

Y- Y +cr Y 

IF ( I T-( I T/2)*2 ) 21,11,21 

11 I F ( JT- ( JT/2)*2) 57,12,57 

12 IF ( KT- ( KT/2 ) *2 ) 57,30,57 

21 I F ( JT— ( JT/2 ) * 2 ) 22,57,22 

22 IF(KT-(KT/2)*2) 30,57,30 
30 RX ( M ) = X 

RY(M) =Y+YLAX( J ) 

RZ ( M ) = Z 
M = M+ 1 

IF ( IT.NE. I VAC X ) GO TO 57 
IF ( JT .NE.IVACY) GO TO 57 
IF (KT.NE. IVACZ) GO TO 57 
NVAC=M- 1 

57 IT = I T + l 

58 CONTINUE 
KT=KT+1 

59 CONTINUE 

JT = JT + 1 

I F ( IDEEP-JT ) 60,110, 60 

60 CONTINUE 
L L =M- 1 

100 RETURN 
110 LD=M-1 

GO TO 60 
END 
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;$ojca$e ##3* a*#######* a*:*# 



c 



£ ajc 5$c £ rfc^;^C^5j«5i|c5?C55:5?<5!S5}:^>!C3te>>C^Oi«:5^ >£•%. akafc# 

'POTENTIAL WELL CONFIGURATION' 

# # # # if. Jf if jjc jjt * # :$: if # A if »£ 5}: * £ # * sjc it * if # # * 



SUBROUTINE PLACE OFFSETS THE INTERSTITIAL IN THE 
'RELAXED' CRYSTAL. 



SUBROUTINE PLACE 

COMMON/ COM 1/RX( 500) ,RY(500) , RZ( 500) ,LCUT( 500) , 

1 LL , LD , I TYPE. NV AC 

COMMON /C0M3/RXI { 500) , RY I ( 500 ) , R Z I ( 500 ) ,CVR,EVR, 
1 NT, TIME , DT ,DTI , I L A Y , RXK { 500 ) ,RYK( 500) ,RZK{ 500) 
COMMON/ C0M4/ IX, IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y, 
ID 1 Z , I VAC X, I VACY, I VACZ 
COMMON/ C0M9/XN VAC, YNVAC.ZNVAC 





GO TO 


( 10 


,20,30 


,4 


0) , I 


10 


LCUT { 


NVAC) 


= 1 








LCUT ( 


1 ) = 


1 








RX( 1 ) 


= 0. 0 










RY ( 1 ) 


= 0.0 










RZ ( 1 ) 


=0.0 










GO TO 


50 








20 


RX { 1) 


= RX ( 1 


) +D1 X 








RY (1) 


= RY ( 1 


J+D1Y 








RZ ( 1 ) 


= RZ ( 1 


) +D 1 Z 








GO TO 


50 








3D 


LCUT ( 


NVAC) 


= I 








R X ( 1 ) 


= RX 


(NVAC) 








RYU ) 


= RY 


(NVAC) 








RZ { 1 ) 


= RZ 


( NVAC) 








GO TO 


50 








40 


RX { 1 ) 


= R X (1 


) + DL X 








R Y ( 1) 


=R Y { 1 


)+DLY 








RZ ( 1 ) 


= RZ( 1 


) +D1 Z 






50 


CONTI 


NUE 










RX(NVAC) =R 


X( NVAC 


) + 


XNVAC 




RY { NVAC ) = RY ( N\/ AC 


) + 


YNVAC 




RZ (NVAC) =R 


Z (NV AC 


) + 


ZNVAC 



R X I (NVAC)=RX(NVAC) 
RYI { NVAC ) =RY ( NVAC) 
RZI(NVAC)=RZ(NVAC) 
RXK (NVAC) =RX (NVAC) 
RYK (NVAC)=RY(NVAC) 
RZK( NVAC )=RZ (NVAC) 
RXI ( 1 ) =RX ( 1 ) 

RYI { 1 ) =R Y( I) 

RZ I { 1 ) = RZ ( 1 ) 

RXK ( I ) =R X { I ) 

RYK { I ) = R Y ( I ) 
RZK(I)=RZ(1 ) 

RETURN 

END 



& ## # >K jfe 



•OTHER CONFIGURATIONS' 

£ # # * j(c * £ if if if # ^ * # * # # >*c 5}: if if if if if jjc if if if if if if if :£ 



SUBROUTINE PLACE CREATES A VACANCY OR IMPLANTS AN 
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INTERSTITIAL OR SELF INTERSTITIAL IN THE LATTICE 
LATTICE ATOMS ARE ALSO OFFSET AS REQUIRED BY 
INPUT DATA. 



SUBROUTINE PLACE 



COMMON/ COM 1/RX ( 500 ) , R Y { 5 00 ) , RZ ( 500) ,LCUT( 500) , 

1 LL,LD, I TYPE, NVAC 

COMMON /COM3/RXI ( 500) , RYI (500 ) , RZI ( 500) ,CVR,EVR, 
1NT .TIME , DT ,DTI , ILAY,RXK( 500) ,RYK( 500) ,RZK( 500) 
COMM ON/ COM4 /I X, I Y, I Z , SCX , SCY , S C Z , I DE EP , D IX , D 1Y, 
1D1Z, I VAC X. I VACY, I V ACZ 
COMMON/ C0M9/XN VAC, YNVAC.ZNVAC 
GO TO ( 10,20,30,4 0), I TYPE 

10 LCUT(NVAC) = 1 
LCUT ( 1 ) = 1 

RX( 1) =0. 0 
RY { 1 ) =0 . 0 
RZ ( 1 ) =0 .0 
GO TO 50 

20 RX ( 1 ) = RX ( NV AC ) + D 1 X 
RY ( 1 ) = RY ( NV AC) +D1Y 
RZ( 1)=RZ(NVAC) +D1Z 
GO TO 50 

30 LCUT(NVAC) = 1 
RX ( I ) = RX(NVAC) 

RY C 1 ) = RY(NVAC) 

RZ( 1) = RZ ( NVAC) 

GO TO 50 

40 RX ( 1 ) = R X ( NV AC ) + D 1 X 
RY ( 1) =RY( NVAC) + D1Y 
RZ ( 1 ) = R Z ( NV AC ) +D 1Z 

50 CONTINUE 

RX(NVAC) =RX( NV AC ) + X NVAC 
RY (NVAC ) =RY ( NV AC ) + YNVAC 
RZ(NVAC) =RZ(NVAC)+ZNVAC 
R XI (NVAC)=RX(NVAC) 

RY I ( NVAC ) = RY (NVAC) 

RZI (NVAC) =RZ (NVAC) 

RXK ( NVAC ) = RX ( N VAC ) 

RYK (NVAC )=RY (NVAC) 

RZK(NVAC)=RZ(NVAC) 

RXI ( 1 ) =RX( 1 ) 

RYI (1 )=RY(1 ) 

RZI ( 1) =RZ( 1) 

RXK ( 1 ) = RX ( 1) 

RYK( 1) = RY (1 ) 

R ZK ( 1 ) = RZ ( 1) 

RETURN 

END 

& 3^ 5?C ^ ^ 5}C ## ^ ^ 



SUBROUTINE STEP CALCULATES FORCES ON THE INTERSTITIAL 
AND ALL OTHER ATOMS OF THE CRYSTAL. 



SUBROUTINE STEP 
C 



c 

COMMON/C OM 1/RX ( 500 ) , RY ( 500) , RZ { 500) ,LCUT( 500) , 
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1 LL » L D * ITYPE.NV AC 

COMMON/C CM5 /ROE , R0E2, ROEM, EX A, EXB, PEXA, PEXB, FXA, PFXA, 
1IQ,TSAVE,BSAVE 

COMMON/ COM 6/ F X ( 500) , FY( 500) , FZ( 500) , PAC ,PFPTC ,FM 
CCMMCN/C0M8/R0EA, ROEB, ROEC, R0EC2 , CPO , CPI , CP2, CP 3, 
1CF0,CF1,CF2,C3D1,CGD2,CGB1 ,CGB2 ,CGF1 ,CGF2 
IF ( 10-1 ) 100,101, 102 

100 I P =2 

GO TO 200 

101 IP=1 

GO TO 200 

102 1=1 

I P=2 

105 DO 195 J=I P , LL 

I F ( LCUT ( J ) ) 195,110,195 

110 DRX=RX( J ) -RX ( I ) 

IF(DRX) 113,117,117 
113 I F ( DRX + ROE ) 195, 195,120 
117 IFIDRX-ROE) 120,195,195 
120 DR Y= RY ( J ) -R Y ( I ) 

I F ( DRY ) 123,127,127 
123 I F ( DRY+R OE ) 195,195,130 
127 I F ( DRY-ROE ) 130, 195,195 
130 DR Z=RZ ( J ) — R Z ( I ) 

IF(DRZ) 133,137,137 
133 I F( DRZ+ROE ) 195,195,140 
137 IF(DRZ-ROE) 140,195,195 
140 DI ST=DRX*DRX+DRY*DRY+DRZ*DRZ 
IF ( D IST-R0E2 ) 150, 195, 195 
150 DIST=SQRT(DIST) 

160 I F ( D I ST-ROEM ) 162, 162,165 
162 FORCE=EXP(PFXA+PEXB*DIST) 

GO TO 180 
165 DFF=ROE-DI ST 

IF(DFF— 1 .0E-10) 195,195, 167 

167 FORCE =(EXP( PAC+PEXB*DI ST ) -P FPT C ) / DFF 
180 I F ( FM-FC RCE ) 190,190, 195 
190 F OD = FORCE/DI ST 
FA=FOD*DRX 
FX ( J ) = FX ( J ) +FA 
F X ( I ) = F X ( I ) -FA 
FA=FOD*DRY 
FY ( J ) = FY ( J ) +FA 
FY( I )=FY( I ) — F A 
FA=FOD*DRZ 
FZ ( J ) = F Z ( J ) +FA 
FZ ( I ) =FZ ( I ) —FA 
195 CONTINUE 

200 DO 300 I =IP, LD 

I F ( LCUT ( I) ) 300,205,300 
205 IP=I+1 

DO 295 J = 1 P , LL 
I F ( LCUT ( J ) ) 295,210,295 
210 DRX=RX ( J )-RX ( I ) 

IF(DRX) 213,217,217 
213 I F ( DRX + ROEC ) 295, 295,220 
217 I F ( DRX-ROE C ) 220,295, 295 
220 DR Y=R Y( J )-RY( I ) 

I F ( DRY ) 223,227,227 
223 IF (DRY+ROEC) 295,295,230 
227 I F ( DRY-ROEC ) 230,295,295 
230 DR Z = RZ ( J ) — R Z ( I ) 

IF(DRZ) 233,237,237 
233 I F ( DRZ +ROEC ) 295,295, 240 
237 IF (DRZ-ROEC) 240,295,295 
240 DI ST=DR X*DR X+ DR Y*DR Y+DRZ*DRZ 
I F ( D IST-RO EC2 ) 250, 295,295 
250 DI ST = SQRT ( DI ST ) 

IF ( D IST-ROE A ) 260,255,255 
255 IF ( D IS T-RO E B ) 265,270,270 
260 FORCE=EXP(FXA*EXB*DIST) 
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GO TO 280 

265 F0RCE=DIST*(DIST*CF2+CF1)+CF0 
GO TO 280 

270 F0RCE=EXP(CGF1+CGB1*DI ST ) - E XP ( C GF 2+CGB2*D I ST) 
280 F 00= FORC E/ DI ST 
FA=F0D*DRX 
FX ( J ) =FX ( J ) +F A 
FX ( I ) =FX { I ) — FA 
FA =FOD*DRY 
FY ( J ) = FY { J ) +F A 
FY ( I ) =FY ( I ) —FA 
F A =FOD^DRZ 
FZ ( J ) = F Z ( J ) +F A 
FZ ( I ) = F Z ( I ) — F A 
295 CONTINUE 
300 CONTINUE 
RETURN 
END 



SUBROUTINE ENERGY CALCULATES MUTUAL POTENTIAL ENERGIES 
FOR THE INTERSTITIAL AND ALL OTHER ATOMS OF THE 
CRYSTAL . 

i 



SUBROUTINE ENERGY 



CCMMON/COM 1/RX ( 500 ) , RY( 500) , RZ ( 500) ,LCUT{ 500) , 

1 LL , LD , I TY P c » NY AC 

COMM ON/C 0M5 /ROE , R0E2 » R OE M, EX A , EX B , P EX A , P EX B , FX A , P FX A , 
1 I Q » T SAV E , BS AVE 

COMMON/ C0M7/PPTC,T POT, PP E ( 1000 ) ,T LP E , ROEL , R0EL2, NEW 
C0MMCN/CCM8/R0EA ,ROEB,ROEC , ROEC 2 , C PO , C P 1 ,CP2 ,CP3, 

I CF0,CF1,CF2,CGD1 ,CGD2,CGB1 , C GB 2 , C GF 1 , CGF 2 
IF { I 0-1 ) 100,101,102 

100 I P = 2 

GO TO 200 

101 I P = 1 

GO TO 200 

102 1=1 
IP =2 

105 DO 595 J = I P , LL 

I F { LCUT ( J) ) 595,510,595 
510 DR X=RX( J ) — R X ( I ) 

IF(DRX) 513,517,517 
513 I F ( DRX+ ROE ) 595,595,520 
517 I F ( DRX-ROE ) 520, 595,595 

520 DRY = RY ( J )-RY ( I ) 

I F { DRY ) 523,527,527 

523 I F ( DRY + ROE ) 595,595,530 
527 IF (DRY-ROE) 530,595,595 
530 DRZ=RZ ( J ) -RZ ( I ) 

IF(DRZ) 533,537,537 
533 I F (DRZ+ROE) 595,595,540 
537 I F { DRZ-ROE ) 540, 595,595 
540 DIST=DRX*DRX+DRY*DRY+DRZ*DRZ 
I F ( D I ST-ROE 2 ) 550,595,595 

550 D I ST = SORT ( D I ST ) 

560 POT=EXP(PEXA+PEXB*DIST J-PPTC 
580 TPOT=TPOT+POT 

PPE( I )=PPE( I ) + BS AV E*PO T 
PPE ( J ) =PPE( J) +TSAVE*POT 
595 CONTINUE 
600 CONTINUE 

200 DO 3 00 I =1 P , LD 
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I F ( L CUT ( I ) ) 300,205,300 
205 IP=I+1 

DO 295 J = 1 P , LL 
I F ( LCUT ( J ) ) 295,210,295 
210 DRX = RX( J ) — RX C I ) 

IF(DRX) 213,217,217 
213 I F { DRX +ROE C ) 295,295,220 
217 IF(DRX-ROEC) 220,295,295 
220 DRY = RY ( J ) — RY ( I J 

IF (DRY ) 223,227,227 
223 I F (DRY+ROEC ) 295,295,230 
227 IF(DRY-ROEC) 230,295,295 
230 DRZ=RZ ( J ) — RZ ( I ) 

IF(DRZ) 2 33,237,237 
233 IF ( DRZ + ROEC ) 295,295,240 
237 I F ( DRZ-ROEC ) 240,295,295 
240 DIST=DRX*DRX+DRY*DRY+DRZ*DRZ 
I F ( D I S T — RO E C2 ) 250,295,295 
250 DI ST = SQRT ( DI ST) 

IF C D IST-ROEA) 260,255, 255 
255 IF (DIST-ROEB) 265,270,270 
260 POT= E XP ( E XA+ EX B*D I ST ) 

GO TO 280 

265 POT = DI ST* { D I ST* ( DI S T* C P3 + CP2 ) +CP 1 ) +CPO 
GO TO 280 

2 70 POT=EXP ( CGD1+CGB1*DIST J-EXP I CGD2+CGB2*D I ST) 
280 TPOT=TPOT+POT 
SAVE=0. 5*P0T 
PP E ( I ) = P P E ( I )+SAVE 
PP E ( J) =P PE ( JJ + SAVE 
295 CONTINUE 
300 CONTINUE 
RETURN 
END 



SUBROUTINE PRINT PRINTS PERTINENT GENERAL DATA FOR 
EACH TIMESTEP PRINTOUT. 



C 

C 



SUBROUTINE PRINT 

COMMON/ COM 1 /RX ( 500) ,RY(500) , RZ { 500) ,LCUT( 500) , 

1 LL , LD , ITYPE.NVAC 

COMMON /C 0M2/IH 1 ( 20) ,IH2(8) , IHS ( 10 ) , IHB(6 ) , IHT (6 ) , 
1TARGET ( 4 ) , TMAS , BULLET! 4) , BMAS, PLANE , TEMP , THERM, DDTI I, 
1DDTI F 

COMMON /COM 3/R XI ( 50 0) , RY I ( 5 00 ) , R Z I ( 500 ) , C V R , E VR , 

1NT , T IME , DT , DT I , I L AY , RXK ( 500 ) , RYK ( 500 ) , R ZK ( 500) 

COMM 0N/C0M4/I X, IY, I Z , SCX , SCY ,S CZ , IDEEP,D1X,D1Y, 

1D1Z, I V AC X , I VAC Y , I VACZ 

COMMON/ C0M5/ ROE, R0E2,R0 EM, EX A, E XB , PE XA, P E XB , F XA , P FXA , 
110 ,TSAVE ,BSAVE 

COMMON /C0M8/RCP A ,ROEB,RQEC ,R0EC2,CP0 ,CP 1 , C P2 , CP3 , 
1CF0,CF1,CF2,CGD1 ,CGD2, CGB1, CGB2,CGF1,CGF2 

9710 FORM AT (4 OX, 1 OA 4 , / , 2 8X , 2 0A4 ,/) 

9720 FORMAT ( 9H TARGET - , 4A4 , 10HPR IMARY - , 4A 4, IX, 14H LATTICE 
1 UNIT =,F7.4,4H ANG) 

9730 FORMAT! 4X, 6HMASS = , F7. 2 , 1 3X , 6H MA S S = , F7 . 2 , 9X , 14HL ATT I C 
IE TEMP =F5 . 2 , 7 H DEG K,,18H THERMAL CUTOFF =,F5.2,3H E 
IV/) 

9740 FO RMAT ( 2H I,A4,8H) PLANE, , 18H PRIMARY ENERGY =, 

1 F5.2,21HKEV, CRYSTAL SIZE ( ,12, 3H X ,I2,3H X ,I2,3H 
1 ),, 4X , 16HVAC A NCY IN SITE , 14/) 

9741 FORMAT ( 2H (,A4,8H) PLANE, , 18H PRIMARY ENERGY =, 

1 F5.3,21HKEV, CRYSTAL SIZE ( ,I2,3H X ,I2,3H X 



12 . 
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9742 

9750 

9760 

9765 

9770 

9780 

9791 

401 

402 

403 



13H ),, 4X » I 3 , ' LAYERS ARE FREE TO MOVE',/) 
FORMAT (2H (,A4,8H) PLANE, , 18H PRIMARY ENERGY 
1 F5.2,21HKEV, CRYSTAL SIZE ( ,I2,3H X ,I2,3H 

1 ),, 4X , 2 0HRE P LAC EMENT IN SITE , 14/) 



X , I 2, 3H 



1 

1 

1 

1 



IMPLANTATION 

•Z = * , 12, 5X, • 10=' 
• LATTICE 



AT 

t 



SITE # - , 

, I 2 , 4X 



• , 13, 3X, • X=' 



12, 4X, 

' DDT I I = * ,F6.4,4X, 
5X, 

STI TIAL 

t 



ATOM START POINT' 



, 6A4,3X,5HPEXA=,F9.5,2X,5HPEXB=, 

5,2X,5HFX 



FORMAT ( 

' Y = ' ,12 , 4X , 

' DDT I F= ' ♦ F 6 . 4, //, 

• X = ' ,F5 .2,3X, • Y= • , F5.2,3X, • Z= • , F5.2,6X, ' INTER 
START POINT' ,5X,'X=' ,F 5.2,3 X,' Y = ' ,F5.2,3X,'Z = 
1F5.2, , // ) 

FORMAT ( 12H POTENTIAL 
1 F9 . 5 , 2 X , 5HPFXA = ,F9. 5) 

FORMAT ( 12X, 6A4, 3X, 5HEXA =,F9.5,2X, 5HE X9 = ,F9. 
1A =,F9.5/) 

FORMAT ( ' WH EN ' , F 8. 4 , ' < R < ' 

1TIAL PARAMETERS ARE',//,' 

IF 10. 3, ' , C p 2 = ' , FI 0 . 3 , ' , CP3 
1 El 0. 3, ' , CF1 = ' , E 1 0 . 3 , ' , CF 2 
FORMAT ( ' CUT-OFF AT',F5.2, ', 

1SE POTENTIAL PARAMETERS ARE' 

1 F8 . 4 , ' , CGD2 = ' , F8. 4, ' , CGB1 
1', CGF 1 =' , F 8. 4 , ' , CGF2 =',F8.4,//) 

FORMAT ( 36X , 'DDTI I = ' , F 6 . 4/ ,36 X , • DDT I F = • ,F6.4) 
WRITE ( 6,9710) IHS, IH1 
WRITE ( 6 ,9720) T ARG ET , BU L L ET , CV R 
WRITE ( 6,9730) TM A S , BMA S , TE MP , THERM 
GO TO (401,402,403,402), ITYPE 
WRITE ( 6,9740) P LANE , E VR , I X , I Y , I Z , NVAC 
GO TO 405 



, F 8 . 4 , ' THE MATCH 
CPO =• , F10.3, 
= ' ,F10.3,/» ' 

=• , El 0.3,//) 

WHEN R > ' , F6.3, 

, 8A4,// , 10X , ' 
=',F8.4,', CGB2 



ING POTEN 
', CPI =• 
CFO =• 

' LU, MOR 
CGD1 =', 

= ' , F 8 . 4 , 



WRITE 
GO TO 
WRITE 
WRI TE 
WRITE 
WRITE 
WRITE 



(6,9741 ) 
405 



PLANE, EVR, IX, IY , IZ, I LAY 



6, 9742 ) 
6,976 0) 
6, 9765) 
6,9770) 
6,978 0) 

1CGF1 , CGF2 

WRITE ( 6,9790) 
RETURN 
END 



PLANE,EVR,IX,IY,IZ»NVAC 

IHB,PEXA,PEXB,PFXA 

IHT , E XA , E XB ,FXA 

ROEA, R0EB,CP0,CP1,CP2,CP3,CF 

ROEC , ROEB , I H2 , CGD1 , CGD2 , CGB1 

NT, DTI, TIME, DT 



0,CF1,CF2 
, CGB2 , 






•PLUCK CONFIGURATION' 



SUBROUTINE PLUCK CHOOSES THE ATOMS WHICH WILL 
MAKE UP THE SMALLER CRYSTAL USED IN INITIAL DYNAMIC 
CALCULATIONS. PERTINENT DATA FOR THIS SMALLER 
CRYSTAL AND POSITIONS OF ALL ATOMS OF THIS CRYSTAL 
CAN THEN BE PUNCHED CUT ON DATA CARDS TO BE USED AS 
INPUT TO THE DYNAMIC PROGRAM. 



SUBROUTINE PLUCK 

COMM ON /COM 1 /RX ( 5 00) ,RY(500) ,RZ( 500) ,LCUT(5 00) , 
1LL.LD, ITYPE, NVAC 

C0MM0N/C0M4/IX,IY,I Z , SCX , SCY , S C Z , I DE EP , D1 X , D1Y , D1 Z , 
1 IVACX, I VACY, IVACZ 
C0MM0N/C0M10/I XNEW , IYNEW , I ZN EW , I I 

COMMON/COM11/RXNEWI (250) , R YNEW I ( 2 50 ) , RZNEW I ( 250 ) , 
1KEEP( 250 ) , NNUM ( 250) 

1500 I X N E W = 7 

IYNEW= I VACY+3 
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I ZNEW=7 
NM=5 

NI = 8 

II = 2 
NM=0 
NX= 16 
NI 13=10 
NI 14 = 5 

I F ( IYNEW.E0.3) GO TO 1514 
I F ( IYNEW.EQ.5) GO TO 1514 
1535 DO 1539 I = I I » NM 
NNUM { I ) =NI 
1509 N I = N I + 1 
NI =NI+1 
I 1= I 1+4 
NM=NM+4 

IF (II. LE. NX) GO TO 1505 
NX=NX+9 
NI =NI + NI 14 
MM =MM+ 1 

I F ( IYNEW.EQ.MM) GO TO 1600 

NM=NM-1 

GO TO 1515 

1514 NX=9 
NM=4 
NI 13=4 
NI 14=11 

1515 DO 1520 I = I I , NM 
NNUM ( I ) =NI 

1520 N I =N I + 1 
N I =N I + 2 
I 1=11+3 
NM=KM+3 

I F ( II. LE . NX) GO TO 1515 
NX=NX+ 1 6 
NI =N I+NI 13 
MM=MM+ 1 

IF( IYNEW.EO.MM ) GO TO 1600 
NM=NM+ 1 
GO TO 1505 
1600 11=11-1 

RXNEWI < 1 ) =RX ( 1 ) 

RYNEWH 1 ) = RY ( 1 ) 

RZNEWI (1 ) = RZ ( 1 ) 

KEEP ( 1) =1 
NNUM ( 1 )=1 

1700 DO 1750 1=2,11 

RXNEWI ( I ) =RX( NNUM{ I ) ) 
RYNEWI ( I ) = RY ( NNU M ( I ) ) 
RZNEWI ( I ) =RZ ( NNUM ( I ) ) 

KEEP ( I ) =NNUM ( I ) 

1750 CONTINUE 
RETURN 
END 



BLOCK DATA 

C0MM0N/CCM1/RX( 1000) , RY(IOOO) , RZ { 1000 ), LCUT ( 1000 ) , 
ILL »LD» ITYPEtNYAC 

DATA RX/ 100 0*0 . 0 / , RY/ 1 000*0 .0/ , RZ/ 1033*0 .0/ , 
1LCUT/1000*0. 0/ 

COMMON/ C0M3/RX I ( 1 0 00) , R Y I ( 1 000 ) , R Z I ( 1 000 ) , C VR ,E VR , 

1 NT , T I ME , DT , DT I , I L AY , RXK ( 13 3 3 ) ,RYK( 1303) , RZK( 1303) 
DATA RXI / 1000*0. O/.RYI / 1000*0. 0/ ,RZ I / 100 0*0. 0/ 
COMMON/ CCM6/FX( 1000), FY ( 1000 ),FZ( 1000) , P A C , P FP TC , FM 
DATA FX/ 103 0*0. 0/,FY/ 1 033* 3.0/, FZ/ 1300* 0.0/ 

END 
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• COMPUTER PROGRAM 
FOR DYNAMIC SIMULATIONS 



THE DYNAMIC PROGRAM IS USED TO DETERMINE THE MINIMUM 
ENERGY REQUIRED IN A SPECIFIC DIRECTION TO CAUSE AN INTER- 
STITIAL ATOM TO EXIT FROM THE CRYSTAL. THE COMPUTATIONS 
REQUIRED ARE ESSENTIALLY THE SAME AS THOSE IN THE STATIC 
SIMULATIONS. CONSEOJEN TLY , COMMENTS WILL ONLY BE INCLUDED 
TO POINT OJT SIGNIFICANT DIFFERENCES FROM THE STATIC 
PROGRAM. TWO CONFIGURATIONS OF THE DYNAMIC PROGRAM EXIST: 

(1) THE ENTIRE CRYSTAL CONFIGURATION - THE ENTIRE 250 
ATOM CRYSTAL IS USED IN ALL CALCULATIONS. 

(2) THE PLUCK CONFIGURATION - A SMALLER CRYSTAL DETER- 
MINED BY SUBROUTINE PLUCK IS USED FOR ALL CALCULATIONS. 



D I ME NS I ON VX (1000) , VY ( 1000 ) ,VZ ( 1000 ), PKE( 1000 ) 
DIMENSION D X ( 1000) ,DY( 1000) ,DZ( 1000) , PTE (1000) 

DI MENS ION RXK ( 1000) , R YK ( 1000 ) ,RZK( 1000) 

DIMENSION KEEP (500 ) 

COMMON /COM 1 /R X ( 1 000 ) , R Y ( 1 000 ) , R Z ( 1 000 ) , LC UT ( 1 000 ) , 
1LL,LD, I TYP E , NY AC 

COMMON /C0M2 /I H 1( 20) ,IH2 (8) , IHS (10) , IHB(6 ) , I HT (6 ) , 
1TARGET ( 4) tTMAS , BULLETl 4) , BMAS, PLANE , TEMP , THERM 
COMM ON / COM3 / RX I (1000) , RYI (1000 ) , RZ I ( 1000) , CVR, EVR, 

1NT, TIME ,DT, DTI , I LAY 

COMMON/ COM 4/ I X , IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z 
COMMON/COM5/ROE.ROE2, ROEM, EX A , EX B , P EX A , P EX B, FX A , P FX A , 
1I0,TSAVE»B SAVE tCX.CY.CZ 

COMMON/ C0M6/FX { 1000 ) . FY( 1000 ) , FZ( 1000) , PAC ,P FP TC , FM 
COMMON /COM 7/P PTC ,TPOT, PPE( 1000) ,TLPE,ROEL, R0EL2 . NEW 
COMMON/ C0M8/R0E A, ROEB, ROEC , R0EC2 ,CPO ,CP 1 ,CP2 ,CP3 , 
1CF0.CF 1 , CF2 .CGD1 ,CGC2, CGB1 , CGB2 , CGF1 , CGF2 
COMM ON/C 0M9/RX SA VE , RYS AVE . RZSAVE 
COMMON/COMA/ A(4,5),MCR0 
C 

90 10 FORMAT ( 20A4) 

9020 FORMAT (8A4, 3F8 . 5, 2F5.2 ) 

9030 FORMAT (4A4.3F3.5 ,6 A4 , F6 .2) 

9040 FORMAT ( F6. 5 , 5X , I 5 , 6 1 4 , 3F 5. 2 ) 

90 50 FORMAT (10A4,A4,4I3,F8.4, I4.F5.3) 

C 

9610 FORMAT ( 1H1 ) 

9620 FORMAT (47X, ' SJMMARY OF ATOMS '//, 35X , 8A4, ' , NT ='I4,//, 
1 3 ( ' ATOM POSITION BIND ENERGY '),//) 

9630 F0RMAT(3( I 5 , 3F 6 . 2 , F 8 . 4 , 8X ) ) 

9640 FORMAT (/4X.F10 .3.25H EV . TOTAL KINETIC ENERGY, , F 10. 3, 

127H EV, TOTAL POTENTIAL E NE RGY , F 1 0 . 3 , • EV , REDUCTION',/ 
1/60X, 'RADIUS = ' , F5.2, ) 

96 50 FORMAT ( 1 05X,4HPAGE, 13, /, 1H1 ) 

9660 FORMAT (/ ' ATOM DX DY DZ 

1VX VY VZ KE PE TE',/) 

9670 FORMAT! 1 1 8 , 3F 1 0 . 3 , 3 F 1 0 . 1 , 3 F 10 . 4 ) 

9680 FORMAT ( ' SHARP DT DECREASE' , 2 E 1 0 . 3 ) 



C ’ 'ENTIRE CRYSTAL CONFIGURATION' 



9690 FORMAT ( 14, 3F5. 2, 14) 

9691 FORMAT (9F8. 4) 
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c 



•PLUCK CONFIGURATION' 

*******£* * * * ** ** v:** *** ** ** ****** ** * 



9690 FORMAT ( 14, 3F5.2,4( 14) ) 

9691 FORMAT! I 5»3(1X ,F8.4) ,1X,I5) 

* *********************************************************** 



c 



c 



c 



c 



c 



c 



c 



RUNTM=4*60-20 
START = 0. 01*1 TI ME (XX) 
DO 2 1=1,1000 
RXK ( I ) =0 .0 
RYK( I ) = 0. 0 
RZK ( I ) = 0 .0 
VX ( I )=0 .0 
V Y( I ) =0. 0 
VZ( I ) = 0 .0 
PKE ( I ) =0.0 
PPE ( I )=0.0 
PT E ( I ) = 0 .0 
2 RZI ( I ) =0.0 
I SHUT= 1 
NRUN=0 



READ ( 5,9010) 

READ ( 5,9020 ) 
READ ( 5,9030) 

READ ( 5,9030) 



I H 1 

IH2,DC0N, ALPHA, RE, ROEC, ROEL 
BULLET ,BMAS,PEXA,PEXB, I HB, THERM 
TARGET , TMA S , EXA , EXB , IHT ,TEMP 



ROE 2 =3 . 0 

READ ( 5,9050) I H S , PL AN E , L S , I X , I Y , I Z , C VR , MCRO ,DTI 
ROE=SORT (R0E2) 

ROEM = ROE-DTI 
R0EL2=R0EL*RGEL 
CVE = 1 . 6 OE-19 
C VM= 1. 672E-27 
FM = 1 .0 E-10 
FM2=FM* FM 
CVD=CVR*1.0F-10 
CVED=CVE/CVD 
PTMAS=TMAS*CVM 
PBMAS=BMAS“CVM 
HT MAS=0 .5*PTMAS/CVE 
HBMAS=0. 5*PBMAS/CVE 
VFAC=1 .0 

TS AVE=BMAS/ ( BMAS+TMAS ) 

BSA VE=TMAS/( BMAS+TMAS) 



FXA=ALOG(-EXB*CV ED )+EXA 
PFXA=ALOG( -PEXB^CVED) +PEXA 
PPTC = EXP (PEXA+-PEX3*R0E ) 
PAC=ALCG (CVED) +PEXA 
PFPTC=E XP( PAC«-PEXB*ROE ) 



CGD1=AL0G( DCOM) +2 . 3* ALPHA* RE 
CGD2 = AL OG ( 2. 0*DC0N ) +AL PHA* RE 
CGB1 =-2 .0*ALPHA*CVR 
CGB2=-ALPHA*C\/R 
CGF1=AL0G(-CGB1*CVED) +CGD1 
CGF2=AL0G(-CGB2*CVED)+CGD2 

R0EA=1 . 50/CVR 
ROE B=2 . 0 / CVR 
R0EC2=R0EC*R0EC 



A ( 1 , 1 )=1 .0 

A( 1 ,2) =ROE A 
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A ( 1,3)=R0EA*R3EA 

A(1 ,4) =R0EA**3 

A ( 1,5) =EXP( EXA*EXB*ROEA) 

A(2, 1)=1.0 
A(2,2) = R 0 E 3 
A( 2,3) =ROEB*ROEB 
A( 2 , 4 ) = ROE 3~*3 

A ( 2 , 5 ) = EXP (CGD1+CGB1*R0EB)-EXP ( CGD2+CGB2*R0EB ) 

A( 3, 1) =0. 0 

A( 3, 2 ) =- 1 • 0 

A ( 3 , 3 ) =-2. D* R0 EA 

A( 3,4) =-3. 0*R0EA*R0EA 

A ( 3 , 5 ) = EXP ( FXA+EX3*R0EA)/CVED 

A (4,1) =0.0 

A ( 4, 2) =- 1 . 0 

A(4,3)=-2.0*R3EB 

A ( 4 » 4) =- 3. 3~R0EB#R0EB 

A ( 4,5) = ( EXP(C3F1+CGB1*R0EB)-EXP(CGF2+CGB2*R0EB) ) /CVED 
CALL CROSYM 
CP 0=A ( 1 ,5) 

CP 1=A (2,5) 

C P2 = A ( 3 ♦ 5 ) 

CP3=A( 4,5) 

CF0=-CP I-^CVFD 
CF1=-2.0*CP2*CVED 
CF2=-3. 0«CP3*C VED 

5 READ ( 5,9040 ) EV R , NTT , NS , N D , I P , I DEE P , I T YPE , N VAC , D 1 X , 
1 D 1 Y , D1 Z 

IF(NTT.EO.O) GO TO 9999 
10= ITYPE-1 
NTTS=NTT 
EV=EVR A 1 . 0E+3 



>{c £ £ * 4 : * fc* * *** 5(5 ****:**:(< 5 {:* 5 fc jfcsjt*##*****# 

C 'PLUCK CONFIGURATION' 

tfofc # ric & ^ ^ # 5£ a* # ^ 3$c £ & # ^ 



DATA FOR THE PLUCK CRYSTAL IS READ INTO THE COMPUTER 
AND APPROPRIATE PARAMETERS ARE CORRELATED. THIS 
CREATES THE CRYSTAL TO BE USED IN THE SIMULATION. 



READ (5,9690) L L , D 1 X , D 1 Y, D 1Z , N VAC , I XNE W , I YNE W , I Z NE W 

I I =60 

LL = I I 

IX= IXNEW 

I Y = I YNE VJ 

I Z = I Z N E W 

DO 15 1=1, LL 

READ (5,9691) I , RX ( I ) , RY ( I ) , RZ ( I ) , KEEP ( I ) 

15 CONTINUE 



•ENTIRE CRYSTAL CONFIGURATION' 



DATA FOR THE ENTIRE CRYSTAL IS READ INTO THE COMPUTER. 
THIS CREATES THE CRYSTAL TO BE USED IN THE SIMULATION. 

READ (5,9690) LL , D1X , D 1Y , D 1Z , NV AC 
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DO 15 1=1, LL, 3 
K= I + 1 
J= 1+2 

READ (5,9691) RX ( I ) , RY ( I ) , RZ ( I ) , RX ( K ) , RY ( K ) , R Z ( K ) , 
1 RX ( J ) , RY( J ) ,RZ ( J ) 



30 ILAY=ID5EP 

IF(IDEEP) 35,35,40 
35 LD = LL 
I L AY= I Y 

40 DO 45 I =1 , LL 
R XK ( I ) =R X ( I ) 

RY K ( I) = RY{ I ) 

R Z K ( I ) =RZ ( I ) 

RX I ( I ) = RX ( I ) 

RY I ( I ) = RY ( I ) 

45 RZ I ( I ) =RZ l I ) 

RX I ( 1 ) = RX ( 1 ) 

RY 111)= RY Cl) 

RZI ( 1 ) = RZ ( 1) 

RXK ( 1 ) = RX ( 1 ) 

RYK(l) =RY( 1) 

R ZK ( 1 ) = RZ ( 1) 

C 

TPOT = 0. 0 
CALL ENERGY 
B I ND=-T POT 
TE =TPOT + B I ND 
TPKE=E V 

TE=TP0T+8IND+TPKE 
WRITE ( 6,9610) 

WRITE ( 6,9620) IH2,NT 
DO 70 1=1, LL, 3 
K = I + 1 
J= I + 2 

70 WRITE (6,9630) KEEP ( I ) , RX ( I ) , RY ( I ) , RZ ( I ) , PP E ( I ) , 

1KEEP(K) ,RX(K. )*RY(K) ,PZ(K) ,PPE ( K ) , KE E P ( J ) , RX ( J ) , RY ( J ) , 
1 RZ ( J ) ,PPE( J) 

WRITE ( 6,9640) TPKE , TPOT ,T E, ROEL 
NPAGE= 1 

WRITE ( 6,9650) NP AGE 



APPROPRIATE IMPACT POINTS FOR INTERSTITIAL ESCAPE ARE 
CHOSEN AND AN IMPACT POINT GENERATOR IS CREATED TO 
GENERATE AND THEN TEST POSSIBLE DIRECTIONS OF ESCAPE. 



CX=3 .8 
CY=D. DO 
CZ= 3 . 6 
NTT=NTT S 
DO 3000 II =1 ,4 
CX = CX+0 . 4 
DO 3000 J J =1 , 3 
CZ=CZ+0. 4 
RXS AVE=CX 
RYSAVE=CY 
R ZSA VE =CZ 

AC= SQR T ( (CX-RXI ( 1) )**2+(CY-RYI ( 1) ) **2+ ( CZ-RZ I ( 1) )**2) 
COX= ( CX-RX I (1) ) / AC 
C0Y=CY-RYI (1) ) /AC 
VO L= SORT (EV/HBMAS) 

VX ( 1) =VOL* COX 
VY ( l)=VOL*COY 
COZ= ( CZ— RZ I ( 1 ) ) /AC 
VZ ( 1 ) =VOL*COZ 
IF(NRUN.EO.O) GO TO 60 
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50 DO 55 1=1, LL 
LCUT( I ) =3 
R X ( I ) = RX I ( I ) 

RY ( I ) = RY Id) 

RZ ( I ) =R Z I ( I ) 

RXK ( I ) = RX I ( I ) 

RY K ( I ) = RY I ( I ) 

55 RZK( I) =RZI ( I ) 

60 NRUN = 1 

DO 6 5 I =2 , L L 
VX( I 1 = 0.0 
VY ( I ) = 0 .0 
VZ(I ) =0.0 
PPE ( I) =0. 0 
PKE ( I )=0 .0 
65 PTE( I ) =3.0 
TPOT= 0 . 0 
NEW = 0 

T I ME = 0. 0 
NT = 0 

IF (NRUN.GT.l. 0) GO TO 95 
NPAGE= NP AGE+ 1 

95 TFAC=2. 0*PTMAS*DTI*CVD 
TFACB=2.0*PBMA S*DTI *CVD 
DT=1 .OE-17 
100 DTOD=DT /CVD 

HDT0D=0 . 5*DT03 
DTOM=DT/PTMAS 
HDT0M=0. 5*DTGM 
DTOMB=DT /PBMAS 
HDT0MB=0.5*DT0MB 
200 CALL STEP 

I F ( L CUT ( 1 ) . GT . 0 ) GO TO 240 
1=1 

RXK ( I ) =RX( I ) 

RYK ( I ) = RY ( I ) 

RZK( I) =RZ( I ) 

RX ( I )=RX(I )+DrOD*(HDTOMB*FX( I ) + VX( I ) ) 

RY ( I ) = RY (I )+DTOD*(HDTOMB*FY ( I )+VY( I ) ) 

RZ ( I ) =R Z ( I )+DTOD*(HDTOMB*FZ( I ) +VZ ( I ) ) 

240 DO 245 1=2, LD 

I F { LCUT ( I ) . GT . 0 ) GO TO 245 
RXK ( I ) =RX( I ) 

RYK ( I ) = RY ( I ) 

RZKt 1 ) = R Z ( I ) 

RX ( I ) =R X ( I )+DTOD*(HDTOM*FX( I ) + VX( I ) ) 

RY ( I ) = R Y ( I )+DTQD*(HDTOM*FY( I ) + VY( I ) ) 

RZ( I ) =RZ{I)+DTOD* ( HDTOM*FZ(I)+VZ( I ) ) 

245 CONTINUE 
CALL STEP 
EMA X = 3 . 0 
FMAX = 0 . 0 
TIME=T IME+DT 
NT =NT + 1 

IF(LCUT( 1) .GT.O) GO TO 265 
I =1 

VSS =V X ( I ) 

VX( I )=VSS+HDT3MB*FX( I ) 

RX ( I ) =RXK ( I )+(VX( I ) +VSS )*HDTOD 
VSS= VY ( I ) 

VY ( I ) = VSS+HDT3MB*FY ( I ) 

RY( I ) =RYKd )+( VY< I )+VSS)*HDTOD 
VSS = VZ( I ) 

VZ( I )=VSS+HDTOMB*FZ ( I ) 

RZ( I) =RZK(I )+( VZ(I )+VSS)*HDTOD 

PKE ( I ) =VX( I ) *V X( I ) + VY ( I )*VY( I ) + VZ( I )*VZ( I ) 

E MAX = PKE ( I ) 

260 FX( I ) =0. 0 
FY ( I 1 = 0.0 
F Z ( I ) =0 . 0 
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265 DO 280 1=2, LD 

I F ( LCU7 ( I) .GTiOJGO TO 280 
VSS=VX( I ) 

VX( I ) = VSS+HDTDM=f r FX( I ) 

RX( I ) = RXK ( I ) +( VX ( I )+VSS J^HDTOD 
VS S = VY ( I ) 

V Y ( I ) = VS S+HDTO M*F Y ( I ) 

RY( I) = RYK(I ) + (VY( I )+VSS )*HDT0D 
VSS=VZ( I ) 

VZ( I )=VSS+HDT3M*FZ( I ) 

RZ( I )=RZK( I ) + (VZ ( I ) +VSS )*HDTOD 

PKE( I ) =VX( I )*VX( I > + VY(I)*VY( I > + VZ(I)*VZ( I > 

FX ( I 1 = 0.0 

FY ( I )=0.0 

F Z ( I ) = 0.0 

IF(PKE( I ) .GT.EMAX) EMAX = PKE(I) 

280 CONTINUE 
DTL=DT 

I F ( EMAX . EQ .0 .0 ) GO TO 285 
GO TO 290 
285 DT= 1 . OE- 17 
GO TO 300 

290 DT=OTI*CVD/SQRT( EMAX) 

CT IME=0 ,01*ITI ME ( XX) -START 
300 IF (ISHUT. EO.-I ) GO TO 400 
310 IF(NS-NT) 400,400,100 
400 CALL PRINT 

410 TPOT =0. 0 

DO 450 1=1, LL 
PPE ( I) =0.0 
450 PTE ( I ) =0.0 
CALL ENERGY 
PKE ( 1 ) =HBMAS*PKE(1 ) 

T PKE =PK E ( 1) 

PT E ( 1 ) = PKE ( 1 )+PPE( 1 ) 

DO 62 0 I =2 , LL 

PK E ( I ) =HTMAS*PKE ( I ) 

TPKE=T PK E+PK E( I ) 

620 PTE ( I ) =PKE ( I ) + PPE ( I ) 

TE=TPOT+BIND 
WRITE ( 6,9660) 

DTEST = ( RY( 1J-RYI (1) )**2 
IF (DTEST. GT. 0.01) DTEST= 0.01 
700 DO 750 I =1 , LD 

DX( I )=RX( I )-RX 1(1) 

DY ( I ) = RY ( I ) — Rf 1(1) 

DZ ( I ) =R Z ( I J-RZI ( I ) 

IF ( OX ( I )**2.3E. DTEST) GO TO 720 
IF (DY( I )**2.GE. DTEST) GO TO 720 
IF ( DZ ( I )**2.GE. DTEST) GO TO 720 
IF(PPE ( I ) .GE.-3.0) GO TO 720 
GO TO 750 



C • PLUCK CONFIGURAT ION' 



IN ORDER TO FOLLOW THE LOGIC OF THE DYNAMIC PROGRAM 
THE ATOMS CREATED BY SUBROUTINE PLUCK HAVE BEEN 
RENUMBERED CONSECUTIVELY FROM ATOM NUMBER 1, BUT THE 
ATOM NUMBER OF EACH ATOM AS IT WAS IN THE ORIGINAL 
CRYSTAL HAS BEEN SAVED IN AN ARRAY. THE WRITE COMMANDS 
MUST THEREFORE PRINT THE ARRAY KEEP(I) SO THAT 
PRINTED OUTPUT WILL BE IN A FORM TO ALLOW READY 
COMPARISON WITH THE ORIGINAL CRYSTAL. TO ACCOMPLISH 
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THIS IN THE PLUCK CONFIGURATION, THE STATEMENTS NUM- 
BERED 720 AND 965 SHOULD BE CHANGED TO READ: 



720 WRI TE (6,9670) K EE P ( I ) , DX ( I ) , DY ( I ) , DZ ( I ) , V X (I) , VY ( I ) , 

1 VZ ( I ) , PKE ( I) , PPE ( I ) ,PTE(I ) 

965 WRITE (6,9630) KEEP ( I ) , RX ( I ) , RY ( I ) , RZ ( I ) , PPE ( I ) , 

1KEEP(K) ,RX(K),RY(K) ,RZ(K) , P PE ( K) ,KEEP(J) ,RX(J) ,RY(J) , 

1 RZ ( J ) , PP E( J ) 

# sjt 5(c 5}: % * £ jk it # ^ 5)c * £ £ if. -if. % % 4 * * £ # * ^ * * $ * # £ 5$c 4c £ & * # -if. if -if. * £ if if fc if £ £ * jjt i(t 



720 WRITE ( 6,9670) I , DX ( I ) , DY ( I ) , DZ (I) , VX ( I ) , VY ( I ) , 
1VZ( I ) , PKE ( I ) , PPE ( I ) , PTE ( I ) 

750 CONTINUE 

WRITE ( 6,964-0) T PKE , T POT , T E , ROEL 
WRITE ( 6,9650) NP AGE 
NPAGE = NP AGE+1 
IF(NT-NTT) 790,950,950 
790 N S = N S + N D 
GO TO 100 
950 CONTINUE 
NS = 0 

955 WRITE ( 6,9620) IH2.NT 
DO 965 1=1 , LL, 3 
K = I + 1 
J= I +2 

965 WRITE ( 6,9630 ) I , RX ( I ) , RY ( I ) , R Z ( I ) , PP E ( I ) , K , RX( K ) , 
1RY(K) , R Z ( K ) , P PE ( K) ,J,RX(J) ,RY(J),RZ(J),PPE(J) 

WRITE ( 6,9640) T P K5 , TPOT , T E , ROEL 
WRITE ( 6,9650) NPAGE 



AFTER CALCULATIONS OF THE DYNAMICS ASSOCIATED WITH 
ALL IMPACT POINTS ALONG A (100) DIRECTION, THE IMPACT 
POINT GENERATOR MUST BE INCREMENTED. 



IFUJ.E0.3) CZ = 3 . 6 
1000 IF(ISHUT) 9999,3000,3000 
3000 CONTINUE 
STOP 
END 



SUBROUTINES B100, PLACE, AND PLUCK ARE NOT USED IN 
DYNAMIC SIMULATIONS. SUBROUTINES STEP, CROYSM, AND 
ENERGY ARE USED IN THE SAME FORM AS IN THE STATIC 
SIMULATIONS AND ARE NOT REPEATED HERE. SUBROUTINE 
PRINT APPEARS IN DYNAMIC SIMULATIONS AS FOLLOWS: 



SUBROUTINE PRINT 



C0MM0N/C0M1/RX( 1 00 0 ) , R Y ( 1 000 ) , R Z ( 1 000 ), LCU T ( 1000 ) , 
1LL.LD, ITYPE.NVAC 

COMMON/ COM2 /I HL ( 2 0 ) , I H2 ( 8 ) , I HS ( 10 ) , I HB ( 6 ) , I HT ( 6 ) , 
1TARGET ( 4) , TMA S , B UL L ET ( 4 ) , B MA S , P LA N E , TE M P , T HE RM 
C0MM0N/C0M3/RX I ( 1000 ) , RY I ( 1000 ) , RZ I ( 1000 ) , C VR , E VR , 
1NT, TIME ,DT, DTI , I LAY 
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COMMON /COM 4/ 1 X, IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z 
COMMON/ C OM5 / RO E , RO E2 , ROEM , EX A , EXB, P E XA, P E XB , FXA , P FXA , 

1 10, T SAVE ,3 SAVE ,CX,CY,CZ 

COMM ON /COM 8 /ROE A, ROEB, ROEC , R0EC2 , CPO ,CP 1 ,CP2 ,CP3 , 
lCFOtCFl , CF2.CGD1 ,CGD2, CGB1 , CGB2 , C GFi , CGF 2 
C0MM0N/C0M9/RXSAVE , RYS AVE , RZSAVE 
97 10 FORMAT ( 4 OX, 10A4, /, 2 8X, 2 0A4, / ) 

9720 FORMAT ( 9H TARGET - ,4A4 , 10HPR IMARY - ,4A4 , IX , 14HL ATT I CE 
1 UNIT = , F 7 . 4 , 4 H ANG) 

9730 FORMAT (4X,6HMASS = , F7 . 2 , 13X , 6HMASS = , F 7 . 2 , 9X , 14HL ATT I C 
IE TEMP = F5.2,7H DEG K,,18H THERMAL CUTOFF =,F5.2,3H E 
IV/) 

9740 FORMAT ( 2H (,A4,8H) PLANE, , 18H PRIMARY ENERGY = , 

1 F5.2,21HKEV, CRYSTAL SIZE ( ,I2,3H X ,I2,3H X ,I2,3H 
1 ),, 4X , 16HVAC ANCY IN SITE , 14/) 

9741 FORMAT ( 2H (,A4,8H) PLANE, ,18H PRIMARY ENERGY =, 

1 F6.5,21HKEV, CRYSTAL SIZE ( ,12, 3H X ,I2,3H X ,12, 
13H ),, 4X, 16H INTER STI TI AL ( - , F 5 . 2 , 2H ,- , F 5 . 2 , 2H , + , 

1F5.2,12H) FROM SITE ,14/) 

9742 FORMAT ( 2H (,A4,8H) PLANE,, 18H PRIMARY ENERGY =, 

1 F5.2,21HKEV, CRYSTAL SIZE ( ,12, 3H X ,I2,3H X , 1 2 , 3H 

1 ),, 4X , 2 OHRE PLACEMENT IN SITE , 14/) 

9750 FORMAT (' PRIMARY START POINT (LU) X = ',F5.2,', Y = ' , 
1F5.2,', Z = ' , F5.2, 5X, 13, ' LAYERS ARE FREE TO MOVE’, 
110X, *10 = ' ,12) 

9751 FORMAT (' OFFSET FROM EOUILIB (LU) 0X = ' , F 5 . 2 , ' , OY = • , 
1F5.2,' ,0Z=' , F5. 2, 5X, ' PRIMARY ENERGY IS',F6.3,' KEV',/) 

9760 FORMAT ( 1 2H POTENTIAL , 6A4 , 3 X , 5H PE XA= , F9 . 5 , 2X , 5HP EX B= , 

1 F9 .5,2X,5HPFXA=, F9.5) 

9765 FORM AT ( 12X,6A4,3X,5HEXA = , F9 . 5 , 2X , 5H EX B = , F9 . 5, 2 X , 5 HFX 
1 A =, F9.5/) 

9770 FORMAT ( ' WHEN',F8.4,' <R <',F8.4, ' THE MATCHING POTEN 
1TIAL PARAMETERS ARE',//,' CPO =',F10.3,', CPI =' 

1F10.3,', C P 2 = ' , F 1 0 . 3 , ' * CP 3 =',F10.3,/,' CFO =' 

1E10.3,', C FI =' , El 0.3, ', CF2 = ' , E 1 0 . 3 , / / ) 

9780 FORMAT ( ' CUT-OFF AT',F5.2,', WHEN R > *,F6.3,' LU, MOR 
1SE POTENTIAL PARAMETERS ARE', 8A4,//,10X, • CGD1 =', 

1 F 8 . 4 , ' , CGD2 =' , F8 .4, ' , CGB1 =',F8.4,', CGB2 =',F8.4, 
1', CGF 1 =' , F 8. 4, ' , CGF 2 =',F8.4,//) 

9790 FORMAT (10H TIMESTEP , 1 4 , 2 2X , 6HDT I = , F5.4, 5H LU, 

1 , 22H ELAPSED TIME (SEC) =, E10.4,', NEXT TIMESTEP IS 
1= ' , E10 .4/ ) 

WRITE ( 6,9710) IHS,IH1 
WRITE ( 6 ,9720) T A RGE T , BU L L ET , CVR 
WRITE ( 6,9730) TM AS , BMA S , TE MP , THERM 
GOTO (401,402,403,402), ITYPE 

401 WRITE ( 6,9740) P LANE , E VR , I X , I Y , I Z , NVAC 
GO TO 405 

402 WRITE ( 6,9741 ) P L ANE , E VR , I X , I Y , I Z , D1X , D1 Y , D1 Z, N V AC 
GO TO 405 

403 WRITE ( 6,9742) P L ANE , E VR , I X , I Y , I Z , NVAC 

405 WRITE ( 6,9750) RX I ( 1 ) , RY I ( 1 ) , RZ I ( 1 ) , I L AY , I Q 
WR IT E( 6,9751 ) R XS A V E , R YSA VE , R Z SA VE ,EVR 
WRITE ( 6,9760) I HB, P EX A , P EX B , P FX A 
WRITE ( 6,9765) I HT , E XA , E XB ,F XA 

WRITE ( 6,9770) R OEA , ROEB , CPO ,CP 1 ,CP 2 ,CP 3 ,CF 0 , CF 1 , CF2 
WRITE ( 6,9780) R OEC , RO EB , I H2 , CG D1 , CGD2 , C GB 1 , CGB 2, 
1CGF1,CGF2 

WRITE ( 6,9790) NT , DT I , T I M E , DT 
RETURN 
END 



BLOCK DATA 

COMMON /COM 1/RX( 1000) ,RY( 1000) ,RZ( 1000) ,LCUT(1000) , 

1 LL , LD, ITYPE, NVAC 

DATA RX/ 100 0*0. 0/, RY/ 1 000*0 . 0/ , R Z/ 1 000* 0 . 0 / , 
1LCUT/100 0*0.0/ 

C0MM0N/C0M3/RX I ( 10 00 ) , RY I ( 10 00 ) , R Z I ( 1000 ) ,CVR,EVR, 
1NT, TIME ,DT, DTI , I LAY 

DATA RX 1/1000*0.0/ ,RYI / 1 000* 0 . 0/ , R Z I / 100 0*0. 0/ 
COMMON/ C0M6/FX (1000) , FY ( 1 0 00 ) , FZ ( 1 000 ) , P AC , P FPTC , FM 
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DATA FX/1000*0.0/,FY/1000*0.0/ tFZ/1000*0.0/ 
END 
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